jean perrot / ctools

Commit 7295db57b1058a5bbbcff45625768ce2e545aba8

Authored by chiachun
1 parent 69ab4467

add ctpsfcube

Inline Side-by-side
Showing 12 changed files with 1008 additions and 17 deletions
configure.ac
  View file @7295db5
... ... @@ -883,6 +883,7 @@ AC_CONFIG_FILES([Makefile
883 883 src/ctskymap/Makefile
884 884 src/ctselect/Makefile
885 885 src/ctexpcube/Makefile
  886 + src/ctpsfcube/Makefile
886 887 src/ctcubemask/Makefile
887 888 src/ctools-setup
888 889 test/Makefile
... ...
pyext/ctexpcube.i
  View file @7295db5
... ... @@ -60,7 +60,7 @@ public:
60 60  
61 61  
62 62 /***********************************************************************//**
63   - * @brief CTA exposure map tool Python extension
  63 + * @brief CTA exposure cube tool Python extension
64 64 ***************************************************************************/
65 65 %extend ctexpcube {
66 66 ctexpcube copy() {
... ...
pyext/ctpsfcube.i 0 → 100644
  View file @7295db5
  1 +/***************************************************************************
  2 + * ctpsfcube - CTA PSF cube tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 by Chia-Chun Lu *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctpsfcube.i
  23 + * @brief CTA PSF cube tool definition
  24 + * @author Chia-Chun Lu
  25 + */
  26 +
  27 +%{
  28 +/* Put headers and other declarations here that are needed for compilation */
  29 +#include "ctpsfcube.hpp"
  30 +#include "GTools.hpp"
  31 +%}
  32 +
  33 +/***********************************************************************//**
  34 + * @class ctpsfcube
  35 + *
  36 + * @brief CTA PSF cube tool
  37 + ***************************************************************************/
  38 +class ctpsfcube : public GApplication {
  39 +
  40 +public:
  41 + // Constructors and destructors
  42 + ctpsfcube(void);
  43 + explicit ctpsfcube(GObservations obs);
  44 + ctpsfcube(int argc, char *argv[]);
  45 + ctpsfcube(const ctpsfcube& app);
  46 + virtual ~ctpsfcube(void);
  47 +
  48 + // Operators
  49 + ctpsfcube& operator=(const ctpsfcube& app);
  50 +
  51 + // Methods
  52 + void clear(void);
  53 + void execute(void);
  54 + void run(void);
  55 + void save(void);
  56 + GObservations& obs(void) { return m_obs; }
  57 + void get_parameters(void);
  58 + void get_obs(void);
  59 + void set_response(void);
  60 + void get_ebounds(void);
  61 + void set_from_cntmap(const std::string& filename);
  62 +};
  63 +
  64 +/***********************************************************************//**
  65 + * @brief CTA PSF cube tool Python extension
  66 + ***************************************************************************/
  67 +%extend ctpsfcube {
  68 + ctpsfcube copy() {
  69 + return (*self);
  70 + }
  71 +}
... ...
src/Makefile.am
  View file @7295db5
... ... @@ -12,7 +12,8 @@ SUBDIRS = ctbin \
12 12 ctselect \
13 13 ctskymap \
14 14 ctexpcube \
15   - ctcubemask
  15 + ctpsfcube \
  16 + ctcubemask
16 17  
17 18 # Scripts to be distributed
18 19 scripts = ctools-setup \
... ... @@ -31,6 +32,7 @@ libctools_la_LIBADD = ctbin/libctbin.la \
31 32 ctselect/libctselect.la \
32 33 ctskymap/libctskymap.la \
33 34 ctexpcube/libctexpcube.la \
  35 + ctpsfcube/libctpsfcube.la \
34 36 ctcubemask/libctcubemask.la
35 37  
36 38 # Distribute scripts
... ...
src/ctcubemask/ctcubemask.cpp
  View file @7295db5
... ... @@ -201,7 +201,7 @@ void ctcubemask::clear(void)
201 201 /***********************************************************************//**
202 202 * @brief Execute application
203 203 *
204   - * This method performs the event selection and saves the result
  204 + * This method applies the cube mask and saves the result.
205 205 ***************************************************************************/
206 206 void ctcubemask::execute(void)
207 207 {
... ... @@ -286,7 +286,7 @@ void ctcubemask::run(void)
286 286 // Save event file name (for possible saving)
287 287 m_infiles[i] = obs->eventfile();
288 288  
289   - // Load observation from temporary file, including event selection
  289 + // Apply mask on the event cube
290 290 apply_mask(obs);
291 291  
292 292 } // endif: had a CTA observation
... ... @@ -302,7 +302,7 @@ void ctcubemask::run(void)
302 302 // Write observation(s) into logger
303 303 if (logTerse()) {
304 304 log << std::endl;
305   - log.header1("Observations after selection");
  305 + log.header1("Observations after masked");
306 306 log << m_obs << std::endl;
307 307 }
308 308  
... ... @@ -312,22 +312,22 @@ void ctcubemask::run(void)
312 312  
313 313  
314 314 /***********************************************************************//**
315   - * @brief Save the selected event list(s)
  315 + * @brief Save the masked event cube(s)
316 316 *
317   - * This method saves the selected event list(s) into FITS file(s). There are
  317 + * This method saves the masked event cube(s) into FITS file(s). There are
318 318 * two modes, depending on the m_use_xml flag.
319 319 *
320   - * If m_use_xml is true, all selected event list(s) will be saved into FITS
  320 + * If m_use_xml is true, all masked event cube(s) will be saved into FITS
321 321 * files, where the output filenames are constructued from the input
322 322 * filenames by prepending the m_prefix string to name. Any path information
323   - * will be stripped form the input name, hence event files will be written
  323 + * will be stripped form the input name, hence event cube files will be written
324 324 * into the local working directory (unless some path information is present
325 325 * in the prefix). In addition, an XML file will be created that gathers
326   - * the filename information for the selected event list(s). If an XML file
  326 + * the filename information for the masked event cube(s). If an XML file
327 327 * was present on input, all metadata information will be copied from this
328 328 * input file.
329 329 *
330   - * If m_use_xml is false, the selected event list will be saved into a FITS
  330 + * If m_use_xml is false, the masked event cubes will be saved into a FITS
331 331 * file.
332 332 ***************************************************************************/
333 333 void ctcubemask::save(void)
... ... @@ -362,10 +362,10 @@ void ctcubemask::save(void)
362 362 * @brief Get application parameters
363 363 *
364 364 * Get all task parameters from parameter file or (if required) by querying
365   - * the user. Times are assumed to be in the native CTA MJD format.
  365 + * the user.
366 366 *
367 367 * This method also loads observations if no observations are yet allocated.
368   - * Observations are either loaded from a single CTA even list, or from a
  368 + * Observations are either loaded from a single CTA even cube, or from a
369 369 * XML file using the metadata information that is stored in that file.
370 370 ***************************************************************************/
371 371 void ctcubemask::get_parameters(void)
... ... @@ -383,7 +383,7 @@ void ctcubemask::get_parameters(void)
383 383 // Try first to open as FITS file
384 384 try {
385 385  
386   - // Load event list in CTA observation
  386 + // Load event cube in CTA observation
387 387 obs.load(m_infile);
388 388  
389 389 // Append CTA observation to container
... ... @@ -430,11 +430,13 @@ void ctcubemask::get_parameters(void)
430 430  
431 431  
432 432 /***********************************************************************//**
433   - * @brief Apply mask to counts cube
  433 + * @brief Mask event cube
434 434 *
435 435 * @param[in] obs CTA observation.
436 436 *
437   - * Masks out event bins by setting their content to -1.
  437 + * Mask an event cube. The mask sets pixel values outside the region of interest
  438 + * or the energy of interest to -1. These pixels will be ignored in likelihood
  439 + * fitting.
438 440 ***************************************************************************/
439 441 void ctcubemask::apply_mask(GCTAObservation* obs)
440 442 {
... ...
src/ctexpcube/ctexpcube.cpp
  View file @7295db5
... ... @@ -223,7 +223,7 @@ void ctexpcube::execute(void)
223 223  
224 224  
225 225 /***********************************************************************//**
226   - * @brief Generate the model map(s)
  226 + * @brief Generate the exposure cube(s).
227 227 *
228 228 * This method reads the task parameters from the parfile, sets up the
229 229 * observation container, loops over all CTA observations in the container
... ...
src/ctpsfcube/Makefile.am 0 → 100644
  View file @7295db5
  1 +## Process this file with automake to produce Makefile.in
  2 +
  3 +# Program name
  4 +bin_PROGRAMS = ctpsfcube
  5 +
  6 +# Sources
  7 +sources = ctpsfcube.cpp
  8 +
  9 +# Libtool library
  10 +noinst_LTLIBRARIES = libctpsfcube.la
  11 +libctpsfcube_la_SOURCES = $(sources)
  12 +
  13 +# Program definition
  14 +ctpsfcube_SOURCES = main.cpp
  15 +ctpsfcube_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
  16 +ctpsfcube_CPPFLAGS = @CPPFLAGS@
  17 +ctpsfcube_LDADD = libctpsfcube.la
  18 +
  19 +# Parameter files
  20 +pfilesdir = $(prefix)/syspfiles
  21 +pfiles_DATA = ctpsfcube.par
... ...
src/ctpsfcube/ctpsfcube 0 → 100755
View file @7295db5
No preview for this file type
src/ctpsfcube/ctpsfcube.cpp 0 → 100644
  View file @7295db5
  1 +/***************************************************************************
  2 + * ctpsfcube - CTA PSF cube tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 by Chia-Chun Lu *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctpsfcube.cpp
  23 + * @brief CTA PSF cube tool 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 "ctpsfcube.hpp"
  33 +#include "GTools.hpp"
  34 +#include "GWcs.hpp"
  35 +
  36 +/* __ Method name definitions ____________________________________________ */
  37 +#define G_GET_EBOUNDS "ctpsfcube::get_ebounds()"
  38 +#define G_GET_PARAMETERS "ctpsfcube::get_parameters()"
  39 +#define G_SET_FROM_CNTMAP "ctpsfcube::set_from_cntmap(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 +ctpsfcube::ctpsfcube(void) : GApplication(CTPSFCUBE_NAME, CTPSFCUBE_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 method creates an instance of the class by copying an existing
  72 + * observations container.
  73 + ***************************************************************************/
  74 +ctpsfcube::ctpsfcube(GObservations obs) : GApplication(CTPSFCUBE_NAME, CTPSFCUBE_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 + * @brief Command line constructor
  93 + *
  94 + * @param[in] argc Number of arguments in command line.
  95 + * @param[in] argv Array of command line arguments.
  96 + ***************************************************************************/
  97 +ctpsfcube::ctpsfcube(int argc, char *argv[]) :
  98 + GApplication(CTPSFCUBE_NAME, CTPSFCUBE_VERSION, argc, argv)
  99 +{
  100 + // Initialise members
  101 + init_members();
  102 +
  103 + // Write header into logger
  104 + log_header();
  105 +
  106 + // Return
  107 + return;
  108 +}
  109 +
  110 +
  111 +/***********************************************************************//**
  112 + * @brief Copy constructor
  113 + *
  114 + * @param[in] app ctpsfcube application.
  115 + ***************************************************************************/
  116 +ctpsfcube::ctpsfcube(const ctpsfcube& app) : GApplication(app)
  117 +{
  118 + // Initialise members
  119 + init_members();
  120 +
  121 + // Copy members
  122 + copy_members(app);
  123 +
  124 + // Return
  125 + return;
  126 +}
  127 +
  128 +
  129 +/***********************************************************************//**
  130 + * @brief Destructor
  131 + ***************************************************************************/
  132 +ctpsfcube::~ctpsfcube(void)
  133 +{
  134 + // Free members
  135 + free_members();
  136 +
  137 + // Return
  138 + return;
  139 +}
  140 +
  141 +
  142 +/*==========================================================================
  143 + = =
  144 + = Operators =
  145 + = =
  146 + ==========================================================================*/
  147 +
  148 +/***********************************************************************//**
  149 + * @brief Assignment operator
  150 + *
  151 + * @param[in] app ctpsfcube application.
  152 + * @return Returns ctpsfcube application.
  153 + ***************************************************************************/
  154 +ctpsfcube& ctpsfcube::operator=(const ctpsfcube& 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 ctpsfcube::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 is the main execution method of the ctpsfcube class. It is invoked
  206 + * when the executable is called from command line. The method generates
  207 + * the PSF cube and saves the result.
  208 + ***************************************************************************/
  209 +void ctpsfcube::execute(void)
  210 +{
  211 + // Signal that some parameters should be read ahead
  212 + m_read_ahead = true;
  213 +
  214 + // Create the PSF cube
  215 + run();
  216 +
  217 + // Save the PSF cube into the output FITS file
  218 + save();
  219 +
  220 + // Return
  221 + return;
  222 +}
  223 +
  224 +
  225 +/***********************************************************************//**
  226 + * @brief Generate the model map(s)
  227 + *
  228 + * This method reads the task parameters from the parfile, sets up the
  229 + * observation container, loops over all CTA observations in the container
  230 + * and generates a PSF cube from the CTA observations.
  231 + ***************************************************************************/
  232 +void ctpsfcube::run(void)
  233 +{
  234 + // If we're in debug mode then all output is also dumped on the screen
  235 + if (logDebug()) {
  236 + log.cout(true);
  237 + }
  238 +
  239 + // Get task parameters
  240 + get_parameters();
  241 +
  242 + // Write parameters into logger
  243 + if (logTerse()) {
  244 + log_parameters();
  245 + log << std::endl;
  246 + }
  247 +
  248 + // Set energy dispersion flag for all CTA observations and save old
  249 + // values in save_edisp vector
  250 + std::vector<bool> save_edisp;
  251 + save_edisp.assign(m_obs.size(), false);
  252 + for (int i = 0; i < m_obs.size(); ++i) {
  253 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  254 + if (obs != NULL) {
  255 + save_edisp[i] = obs->response()->apply_edisp();
  256 + obs->response()->apply_edisp(m_apply_edisp);
  257 + }
  258 + }
  259 +
  260 + // Write observation(s) into logger
  261 + if (logTerse()) {
  262 + log << std::endl;
  263 + if (m_obs.size() > 1) {
  264 + log.header1("Observations");
  265 + }
  266 + else {
  267 + log.header1("Observation");
  268 + }
  269 + log << m_obs << std::endl;
  270 + }
  271 +
  272 + // Write header
  273 + if (logTerse()) {
  274 + log << std::endl;
  275 + log.header1("Generate PSF cube");
  276 + }
  277 +
  278 + // Fill PSF
  279 + m_psfcube.fill(m_obs);
  280 +
  281 + // Restore energy dispersion flag for all CTA observations
  282 + for (int i = 0; i < m_obs.size(); ++i) {
  283 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  284 + if (obs != NULL) {
  285 + obs->response()->apply_edisp(save_edisp[i]);
  286 + }
  287 + }
  288 +
  289 + // Return
  290 + return;
  291 +}
  292 +
  293 +
  294 +/***********************************************************************//**
  295 + * @brief Save PSF cube
  296 + *
  297 + ***************************************************************************/
  298 +void ctpsfcube::save(void)
  299 +{
  300 + // Write header
  301 + if (logTerse()) {
  302 + log << std::endl;
  303 + log.header1("Save PSF cube");
  304 + }
  305 +
  306 + // Get output filename
  307 + m_outfile = (*this)["outfile"].filename();
  308 +
  309 + // Save PSF cube
  310 + m_psfcube.save(m_outfile, clobber());
  311 +
  312 + // Return
  313 + return;
  314 +}
  315 +
  316 +
  317 +/***********************************************************************//**
  318 + * @brief Get application parameters
  319 + *
  320 + * Get all task parameters from parameter file or (if required) by querying
  321 + * the user. The parameters are read in the correct order.
  322 + *
  323 + * @todo Setup PSF cube from counts map
  324 + ***************************************************************************/
  325 +void ctpsfcube::get_parameters(void)
  326 +{
  327 + // If we do not have any observations in the container then get an
  328 + // input file name or observation descriptor file
  329 + if (m_obs.size() == 0) {
  330 + get_obs();
  331 + }
  332 +
  333 + // Make sure that response is set
  334 + set_response();
  335 +
  336 + // Read energy dispersion flag
  337 + m_apply_edisp = (*this)["edisp"].boolean();
  338 +
  339 + // If no counts map is specified then setup the PSF cube from
  340 + // the user parameters
  341 + std::string cntmap = (*this)["cntmap"].filename();
  342 + if ((cntmap == "NONE") || (gammalib::strip_whitespace(cntmap) == "")) {
  343 +
  344 + // Get user parameters for counts map definition
  345 + std::string wcs = (*this)["proj"].string();
  346 + std::string coordsys = (*this)["coordsys"].string();
  347 + double xref = (*this)["xref"].real();
  348 + double yref = (*this)["yref"].real();
  349 + double binsz = (*this)["binsz"].real();
  350 + int nxpix = (*this)["nxpix"].integer();
  351 + int nypix = (*this)["nypix"].integer();
  352 + double dmax = (*this)["dmax"].real();
  353 + int ndbins = (*this)["ndbins"].integer();
  354 +
  355 + // Get energy definition
  356 + get_ebounds();
  357 +
  358 + // Define exposure cube
  359 + m_psfcube = GCTAMeanPsf(wcs, coordsys, xref, yref,
  360 + -binsz, binsz, nxpix, nypix,
  361 + m_ebounds, dmax, ndbins);
  362 + }
  363 +
  364 + // ... otherwise setup the exposure cube from the counts map
  365 + else {
  366 +
  367 + // Set PSF cube from counts map
  368 + set_from_cntmap(cntmap);
  369 +
  370 + }
  371 +
  372 + // Read output filename (if needed)
  373 + if (m_read_ahead) {
  374 + m_outfile = (*this)["outfile"].filename();
  375 + }
  376 +
  377 + // Return
  378 + return;
  379 +}
  380 +
  381 +
  382 +/***********************************************************************//**
  383 + * @brief Get observation definition
  384 + *
  385 + * Get observation definition from the user parameters.
  386 + ***************************************************************************/
  387 +void ctpsfcube::get_obs(void)
  388 +{
  389 + // Get input filename
  390 + std::string filename = (*this)["infile"].filename();
  391 +
  392 + // Try first to open as FITS file
  393 + try {
  394 +
  395 + // Allocate CTA observation
  396 + GCTAObservation obs;
  397 +
  398 + // Load input file in CTA observation
  399 + obs.load(filename);
  400 +
  401 + // Append CTA observation to container
  402 + m_obs.append(obs);
  403 +
  404 +
  405 + }
  406 +
  407 + // ... otherwise try to open as XML file
  408 + catch (GException::fits_open_error &e) {
  409 +
  410 + // Load observations from XML file
  411 + m_obs.load(filename);
  412 +
  413 + }
  414 +
  415 + // Return
  416 + return;
  417 +}
  418 +
  419 +
  420 +/***********************************************************************//**
  421 + * @brief Set observation response
  422 + *
  423 + * Set response for all observations that so have no response.
  424 + ***************************************************************************/
  425 +void ctpsfcube::set_response(void)
  426 +{
  427 + // Loop over all observations
  428 + for (int i = 0; i < m_obs.size(); ++i) {
  429 +
  430 + // Is this observation a CTA observation?
  431 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  432 +
  433 + // Yes ...
  434 + if (obs != NULL) {
  435 +
  436 + // Set response if we don't have one
  437 + if (!obs->hasresponse()) {
  438 +
  439 + // Load response information
  440 + std::string database = (*this)["caldb"].string();
  441 + std::string irf = (*this)["irf"].string();
  442 +
  443 + // Set calibration database. If specified parameter is a
  444 + // directory then use this as the pathname to the calibration
  445 + // database. Otherwise interpret this as the instrument name,
  446 + // the mission being "cta"
  447 + GCaldb caldb;
  448 + if (gammalib::dir_exists(database)) {
  449 + caldb.rootdir(database);
  450 + }
  451 + else {
  452 + caldb.open("cta", database);
  453 + }
  454 +
  455 + // Set reponse
  456 + obs->response(irf, caldb);
  457 +
  458 + } // endif: observation already has a response
  459 +
  460 + } // endif: observation was a CTA observation
  461 +
  462 + } // endfor: looped over all observations
  463 +
  464 + // Return
  465 + return;
  466 +}
  467 +
  468 +
  469 +/***********************************************************************//**
  470 + * @brief Get energy boundaries from parameters
  471 + *
  472 + * @exception GException::invalid_value
  473 + * Invalid extension name encountered.
  474 + *
  475 + * Get the energy boundaries from the user parameters.
  476 + ***************************************************************************/
  477 +void ctpsfcube::get_ebounds(void)
  478 +{
  479 + // Determine the energy binning alogrithm
  480 + std::string ebinalg = (*this)["ebinalg"].string();
  481 +
  482 + // If we have the binning given by a file then try to get the boundaries
  483 + // from that file
  484 + if (ebinalg == "FILE") {
  485 +
  486 + // Get filename
  487 + std::string filename = (*this)["ebinfile"].filename();
  488 +
  489 + // Open fits file to check which extension is given
  490 + GFits file(filename);
  491 +
  492 + // Check first for EBOUNDS extension
  493 + if (file.contains("EBOUNDS")) {
  494 + file.close();
  495 + m_ebounds.load(filename, "EBOUNDS");
  496 + }
  497 +
  498 + // ... then check for ENERGYBINS extension
  499 + else if (file.contains("ENERGYBINS")) {
  500 + file.close();
  501 + m_ebounds.load(filename, "ENERGYBINS");
  502 + }
  503 +
  504 + // ... otherwise throw an exception
  505 + else {
  506 + file.close();
  507 + std::string msg = "No extension with name \"EBOUNDS\" or"
  508 + " \"ENERGYBINS\" found in FITS file"
  509 + " \""+filename+"\".\n"
  510 + "An \"EBOUNDS\" or \"ENERGYBINS\" extension"
  511 + " is required if the parameter \"ebinalg\""
  512 + " is set to \"FILE\".";
  513 + throw GException::invalid_value(G_GET_EBOUNDS, msg);
  514 + }
  515 + }
  516 +
  517 + // ... otherwise read emin, emax and nebins
  518 + else {
  519 +
  520 + // Get the relevant parameters
  521 + double emin = (*this)["emin"].real();
  522 + double emax = (*this)["emax"].real();
  523 + int enumbins = (*this)["enumbins"].integer();
  524 + bool log = ((*this)["ebinalg"].string() == "LIN") ? false : true;
  525 +
  526 + // Create energy boundaries
  527 + m_ebounds = GEbounds(enumbins,
  528 + GEnergy(emin, "TeV"),
  529 + GEnergy(emax, "TeV"),
  530 + log);
  531 +
  532 + }
  533 +
  534 + // Return
  535 + return;
  536 +}
  537 +
  538 +
  539 +/***********************************************************************//**
  540 + * @brief Set exposure cube definition from counts map
  541 + *
  542 + * @exception GException::invalid_value
  543 + * Invalid counts map projection or invalid events encountered.
  544 + *
  545 + * Set exposure cube definition from counts map.
  546 + ***************************************************************************/
  547 +void ctpsfcube::set_from_cntmap(const std::string& filename)
  548 +{
  549 + // Allocate CTA observation
  550 + GCTAObservation obs;
  551 +
  552 + // Load counts map in CTA observation
  553 + obs.load(filename);
  554 +
  555 + // Set exposure cube from counts map
  556 + const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(obs.events());
  557 +
  558 + // Continue only if cube is valid
  559 + if (cube != NULL) {
  560 +
  561 + // Get sky map projection
  562 + const GWcs* wcs = dynamic_cast<const GWcs*>(cube->map().projection());
  563 +
  564 + // Continue only if projection is valid
  565 + if (wcs != NULL) {
  566 +
  567 + // Get user parameters for counts map definition
  568 + std::string proj = wcs->code();
  569 + std::string coordsys = wcs->coordsys();
  570 + double xref = wcs->crval(0);
  571 + double yref = wcs->crval(1);
  572 + double dx = wcs->cdelt(0);
  573 + double dy = wcs->cdelt(1);
  574 + int nx = cube->map().nx();
  575 + int ny = cube->map().ny();
  576 + double dmax = (*this)["dmax"].real();
  577 + int ndbins = (*this)["ndbins"].integer();
  578 +
  579 + // Get energy definition
  580 + m_ebounds = cube->ebounds();
  581 +
  582 + // Define exposure cube
  583 + m_psfcube = GCTAMeanPsf(proj, coordsys, xref, yref,
  584 + dx, dy, nx, ny,
  585 + m_ebounds, dmax, ndbins);
  586 +
  587 + } // endif: WCS projection was valid
  588 +
  589 + // ... projection is not of WCS type
  590 + else {
  591 + std::string msg = "Counts map project is not of WCS type.";
  592 + throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
  593 + }
  594 +
  595 + } // endif: observation contained an events cube
  596 +
  597 + // ... there is not events cube
  598 + else {
  599 + std::string msg = "No events cube found in file \""
  600 + ""+filename+"\".";
  601 + throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
  602 + }
  603 +
  604 + // Return
  605 + return;
  606 +}
  607 +
  608 +
  609 +/*==========================================================================
  610 + = =
  611 + = Private methods =
  612 + = =
  613 + ==========================================================================*/
  614 +
  615 +/***********************************************************************//**
  616 + * @brief Initialise class members
  617 + ***************************************************************************/
  618 +void ctpsfcube::init_members(void)
  619 +{
  620 + // Initialise members
  621 + m_outfile.clear();
  622 + m_apply_edisp = false;
  623 +
  624 + // Initialise protected members
  625 + m_read_ahead = false;
  626 + m_obs.clear();
  627 + m_psfcube.clear();
  628 + m_ebounds.clear();
  629 +
  630 + // Set logger properties
  631 + log.date(true);
  632 +
  633 + // Return
  634 + return;
  635 +}
  636 +
  637 +
  638 +/***********************************************************************//**
  639 + * @brief Copy class members
  640 + *
  641 + * @param[in] app Application.
  642 + ***************************************************************************/
  643 +void ctpsfcube::copy_members(const ctpsfcube& app)
  644 +{
  645 + // Copy attributes
  646 + m_outfile = app.m_outfile;
  647 + m_apply_edisp = app.m_apply_edisp;
  648 +
  649 + // Copy protected members
  650 + m_read_ahead = app.m_read_ahead;
  651 + m_obs = app.m_obs;
  652 + m_psfcube = app.m_psfcube;
  653 + m_ebounds = app.m_ebounds;
  654 +
  655 + // Return
  656 + return;
  657 +}
  658 +
  659 +
  660 +/***********************************************************************//**
  661 + * @brief Delete class members
  662 + ***************************************************************************/
  663 +void ctpsfcube::free_members(void)
  664 +{
  665 + // Write separator into logger
  666 + if (logTerse()) {
  667 + log << std::endl;
  668 + }
  669 +
  670 + // Return
  671 + return;
  672 +}
... ...
src/ctpsfcube/ctpsfcube.hpp 0 → 100644
  View file @7295db5
  1 +/***************************************************************************
  2 + * ctpsfcube - CTA PSF cube tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 by Chia-Chun Lu *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctpsfcube.hpp
  23 + * @brief CTA PSF cube tool definition
  24 + * @author Chia-Chun Lu
  25 + */
  26 +
  27 +#ifndef CTPSFCUBE_HPP
  28 +#define CTPSFCUBE_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include "GammaLib.hpp"
  32 +#include "GCTALib.hpp"
  33 +
  34 +/* __Definitions _________________________________________________________ */
  35 +#define CTPSFCUBE_NAME "ctpsfcube"
  36 +#define CTPSFCUBE_VERSION "00-01-00"
  37 +
  38 +
  39 +/***********************************************************************//**
  40 + * @class ctpsfcube
  41 + *
  42 + * @brief CTA PSF cube tool
  43 + ***************************************************************************/
  44 +class ctpsfcube : public GApplication {
  45 +
  46 +public:
  47 + // Constructors and destructors
  48 + ctpsfcube(void);
  49 + explicit ctpsfcube(GObservations obs);
  50 + ctpsfcube(int argc, char *argv[]);
  51 + ctpsfcube(const ctpsfcube& app);
  52 + virtual ~ctpsfcube(void);
  53 +
  54 + // Operators
  55 + ctpsfcube& operator=(const ctpsfcube& app);
  56 +
  57 + // Methods
  58 + void clear(void);
  59 + void execute(void);
  60 + void run(void);
  61 + void save(void);
  62 + GObservations& obs(void) { return m_obs; }
  63 + void get_parameters(void);
  64 + void get_obs(void);
  65 + void set_response(void);
  66 + void get_ebounds(void);
  67 + void set_from_cntmap(const std::string& filename);
  68 +
  69 +protected:
  70 + // Protected methods
  71 + void init_members(void);
  72 + void copy_members(const ctpsfcube& app);
  73 + void free_members(void);
  74 +
  75 + // User parameters
  76 + std::string m_outfile; //!< Output PSF cube file
  77 + bool m_apply_edisp; //!< Apply energy dispersion?
  78 +
  79 + // Protected members
  80 + bool m_read_ahead; //!< Read ahead parameters
  81 + GObservations m_obs; //!< Observation container
  82 + GCTAMeanPsf m_psfcube; //!< PSF cube
  83 + GEbounds m_ebounds; //!< Energy boundaries
  84 +};
  85 +
  86 +#endif /* CTPSFCUBE_HPP */
... ...
src/ctpsfcube/ctpsfcube.par 0 → 100644
  View file @7295db5
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: ctexpcube.par #
  8 +# Version: 00-01-00 #
  9 +# #
  10 +# Author: Juergen Knoedlseder #
  11 +# Jurgen.Knodlseder@irap.omp.eu #
  12 +# IRAP #
  13 +# #
  14 +# Purpose: Parameter file for the ctexpcube application #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# General parameters
  19 +#===================
  20 +infile, f, a, "NONE",,, "Observation definition file"
  21 +cntmap, f, a, "NONE",,, "Counts map for exposure cube definition"
  22 +outfile, f, a, "expcube.fits",,, "Output exposure cube file"
  23 +caldb, s, a, "dummy",,, "Calibration database"
  24 +irf, s, a, "cta_dummy_irf",,, "Instrument response function"
  25 +edisp, b, h, no,,, "Apply energy dispersion?"
  26 +#
  27 +# Energy binning parameters
  28 +#==========================
  29 +ebinalg, s, h, "LOG", FILE|LIN|LOG,,"Algorithm for defining energy bins"
  30 +emin, r, a, 0.1,,,"Start value for first energy bin in TeV"
  31 +emax, r, a, 100.0,,,"Stop value for last energy bin in TeV"
  32 +enumbins, i, a, 20,,,"Number of energy bins"
  33 +ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition"
  34 +#
  35 +# Spatial binning parameters
  36 +#===========================
  37 +nxpix, i, a, 200,,, "Size of the X axis in pixels"
  38 +nypix, i, a, 200,,, "Size of the Y axis in pixels"
  39 +binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
  40 +coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
  41 +xref, r, a, 83.63,,, "First coordinate of image center in degrees (RA or galactic l)"
  42 +yref, r, a, 22.01,,, "Second coordinate of image center in degrees (DEC or galactic b)"
  43 +axisrot, r, h, 0.0,,, "Rotation angle of image axis, in degrees"
  44 +proj, s, a, "CAR", AIT|AZP|CAR|MER|STG|TAN,, "Projection method e.g. AIT|AZP|CAR|MER|STG|TAN"
  45 +#
  46 +# Histogram binning parameters
  47 +#=============================
  48 +dmax, r, a, 0.3,,, "Upper bound of delta"
  49 +ndbins, i, a, 60,,, "Number of delta bins"
  50 +#
  51 +# Standard parameters
  52 +#====================
  53 +chatter, i, h, 2,0,4, "Chattiness of output"
  54 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  55 +debug, b, h, no,,, "Debugging mode activated"
  56 +mode, s, h, "ql",,, "Mode of automatic parameters"
  57 +logfile, f, h, "ctexpcube.log",,, "Log filename"
... ...
src/ctpsfcube/main.cpp 0 → 100644
  View file @7295db5
  1 +/***************************************************************************
  2 + * ctpsfcube - CTA PSF cube tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 by Chia-Chun Lu *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctpsfcube/main.cpp
  23 + * @brief CTA PSF 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 "ctpsfcube.hpp"
  33 +
  34 +
  35 +/***********************************************************************//**
  36 + * @brief Main entry point of application
  37 + *
  38 + * @param[in] argc Number of command line arguments.
  39 + * @param[in] argv Command line arguments.
  40 + *
  41 + * This is the main entry point of the ctpsfcube application. It allocates a
  42 + * ctpsfcube object and runs the application.
  43 + ***************************************************************************/
  44 +int main (int argc, char *argv[])
  45 +{
  46 + // Initialise return code
  47 + int rc = 1;
  48 +
  49 + // Create instance of application
  50 + ctpsfcube application(argc, argv);
  51 +
  52 + // Run application
  53 + try {
  54 + // Execute application
  55 + application.execute();
  56 +
  57 + // Signal success
  58 + rc = 0;
  59 + }
  60 + catch (std::exception &e) {
  61 +
  62 + // Extract error message
  63 + std::string message = e.what();
  64 + std::string signal = "*** ERROR encounterted in the execution of"
  65 + " ctpsfcube. Run aborted ...";
  66 +
  67 + // Write error in logger
  68 + application.log << signal << std::endl;
  69 + application.log << message << std::endl;
  70 +
  71 + // Write error on standard output
  72 + std::cout << signal << std::endl;
  73 + std::cout << message << std::endl;
  74 +
  75 + } // endcatch: catched any application error
  76 +
  77 + // Return
  78 + return rc;
  79 +}
... ...