Veronika Vodeb / ctools

Commit 6f2f518db3dac89ae204497c99a0c9c7e147a828

Authored by Josh Cardenzana
1 parent 18b7c535

Updated ctphase tool (#1982) 

Updated build capabilities (configure.ac and Makefile.am files).
Updated ctphase to add phase information.
Updated GammaLib class "GCTAEventList" to allow specifying whether phase information and detector coordinates are set.
Inline Side-by-side
Showing 7 changed files with 1313 additions and 2 deletions
configure.ac
  View file @6f2f518
... ... @@ -1018,6 +1018,7 @@ AC_CONFIG_FILES([Makefile
1018 1018 src/ctobssim/Makefile
1019 1019 src/ctskymap/Makefile
1020 1020 src/ctselect/Makefile
  1021 + src/ctphase/Makefile
1021 1022 src/ctexpcube/Makefile
1022 1023 src/ctpsfcube/Makefile
1023 1024 src/ctedispcube/Makefile
... ...
src/Makefile.am
  View file @6f2f518
... ... @@ -40,7 +40,8 @@ SUBDIRS = support \
40 40 cttsmap \
41 41 ctbutterfly \
42 42 ctulimit \
43   - cterror
  43 + cterror \
  44 + ctphase
44 45  
45 46 # Scripts to be distributed
46 47 scripts = ctools-setup \
... ... @@ -68,7 +69,9 @@ libctools_la_LIBADD = support/libsupport.la \
68 69 cttsmap/libcttsmap.la \
69 70 ctbutterfly/libctbutterfly.la \
70 71 ctulimit/libctulimit.la \
71   - cterror/libcterror.la
  72 + cterror/libcterror.la \
  73 + ctphase/libctphase.la
  74 +
72 75 libctools_la_LDFLAGS = -version-info $(CTOOLS_LT_VERSION) @OPENMP_LIB_LDFLAGS@
73 76  
74 77 # Distribute scripts
... ...
src/ctphase/Makefile.am 0 → 100644
  View file @6f2f518
  1 +#############################################################################
  2 +# Makefile.am for ctphase tool #
  3 +# ------------------------------------------------------------------------- #
  4 +# Copyright (C) 2017 Leonardo Di Venere #
  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 +# Set include directories
  23 +INCLUDES = -I$(top_srcdir)/src/support
  24 +
  25 +# Program name
  26 +bin_PROGRAMS = ctphase
  27 +
  28 +# Sources
  29 +sources = ctphase.cpp ctphase.hpp
  30 +
  31 +# Libtool library (only needed for python bindings)
  32 +noinst_LTLIBRARIES = libctphase.la
  33 +libctphase_la_SOURCES = $(sources)
  34 +libctphase_la_LDFLAGS = @OPENMP_LIB_LDFLAGS@
  35 +
  36 +# Program definition
  37 +ctphase_SOURCES = main.cpp
  38 +ctphase_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
  39 +ctphase_CPPFLAGS = @CPPFLAGS@
  40 +ctphase_LDADD = libctphase.la $(top_builddir)/src/support/libsupport.la
  41 +
  42 +# Parameter files
  43 +pfilesdir = $(prefix)/syspfiles
  44 +pfiles_DATA = $(srcdir)/ctphase.par
  45 +
  46 +# Extra files to be distributed in the tar file
  47 +EXTRA_DIST = $(pfiles_DATA)
... ...
src/ctphase/ctphase.cpp 0 → 100644
  View file @6f2f518
  1 +/***************************************************************************
  2 + * ctphase - Data selection tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2010-2017 by Juergen Knoedlseder *
  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 ctphase.cpp
  23 + * @brief Data selection tool definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdlib>
  32 +#include "ctphase.hpp"
  33 +#include "GTools.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +#define G_RUN "ctphase::run()"
  37 +#define G_SELECT_EVENTS "ctphase::select_events(GCTAObservation*, "\
  38 + "std::string&, std::string&, std::string&)"
  39 +#define G_SET_EBOUNDS "ctphase::set_ebounds(GCTAObservation*, GEbounds&)"
  40 +
  41 +/* __ Debug definitions __________________________________________________ */
  42 +
  43 +/* __ Coding definitions _________________________________________________ */
  44 +//#define G_USE_MKSTEMP //!< Use mkstemp for temporary filename
  45 +
  46 +
  47 +/*==========================================================================
  48 + = =
  49 + = Constructors/destructors =
  50 + = =
  51 + ==========================================================================*/
  52 +
  53 +/***********************************************************************//**
  54 + * @brief Void constructor
  55 + ***************************************************************************/
  56 +ctphase::ctphase(void) : ctobservation(CTPHASE_NAME, CTPHASE_VERSION)
  57 +{
  58 + // Initialise members
  59 + init_members();
  60 +
  61 + // Return
  62 + return;
  63 +}
  64 +
  65 +
  66 +/***********************************************************************//**
  67 + * @brief Observations constructor
  68 + *
  69 + * param[in] obs Observation container.
  70 + *
  71 + * Creates an instance of the class that is initialised using the information
  72 + * provided in an observation container.
  73 + ***************************************************************************/
  74 +ctphase::ctphase(const GObservations& obs) :
  75 + ctobservation(CTPHASE_NAME, CTPHASE_VERSION, obs)
  76 +{
  77 + // Initialise members
  78 + init_members();
  79 +
  80 + // Return
  81 + return;
  82 +}
  83 +
  84 +
  85 +
  86 +
  87 +/***********************************************************************//**
  88 + * @brief Command line constructor
  89 + *
  90 + * @param[in] argc Number of arguments in command line.
  91 + * @param[in] argv Array of command line arguments.
  92 + ***************************************************************************/
  93 +ctphase::ctphase(int argc, char *argv[]) :
  94 + ctobservation(CTPHASE_NAME, CTPHASE_VERSION, argc, argv)
  95 +{
  96 + // Initialise members
  97 + init_members();
  98 +
  99 + // Return
  100 + return;
  101 +}
  102 +
  103 +
  104 +/***********************************************************************//**
  105 + * @brief Copy constructor
  106 + *
  107 + * @param[in] app Application.
  108 + ***************************************************************************/
  109 +ctphase::ctphase(const ctphase& app) : ctobservation(app)
  110 +{
  111 + // Initialise members
  112 + init_members();
  113 +
  114 + // Copy members
  115 + copy_members(app);
  116 +
  117 + // Return
  118 + return;
  119 +}
  120 +
  121 +
  122 +/***********************************************************************//**
  123 + * @brief Destructor
  124 + ***************************************************************************/
  125 +ctphase::~ctphase(void)
  126 +{
  127 + // Free members
  128 + free_members();
  129 +
  130 + // Return
  131 + return;
  132 +}
  133 +
  134 +
  135 +/*==========================================================================
  136 + = =
  137 + = Operators =
  138 + = =
  139 + ==========================================================================*/
  140 +
  141 +/***********************************************************************//**
  142 + * @brief Assignment operator
  143 + *
  144 + * @param[in] app Application.
  145 + * @return Application.
  146 + ***************************************************************************/
  147 +ctphase& ctphase::operator=(const ctphase& app)
  148 +{
  149 + // Execute only if object is not identical
  150 + if (this != &app) {
  151 +
  152 + // Copy base class members
  153 + this->ctobservation::operator=(app);
  154 +
  155 + // Free members
  156 + free_members();
  157 +
  158 + // Initialise members
  159 + init_members();
  160 +
  161 + // Copy members
  162 + copy_members(app);
  163 +
  164 + } // endif: object was not identical
  165 +
  166 + // Return this object
  167 + return *this;
  168 +}
  169 +
  170 +
  171 +/*==========================================================================
  172 + = =
  173 + = Public methods =
  174 + = =
  175 + ==========================================================================*/
  176 +
  177 +/***********************************************************************//**
  178 + * @brief Clear ctphase tool
  179 + *
  180 + * Clears ctphase tool.
  181 + ***************************************************************************/
  182 +void ctphase::clear(void)
  183 +{
  184 + // Free members
  185 + free_members();
  186 + this->ctobservation::free_members();
  187 + this->ctool::free_members();
  188 +
  189 + // Clear base class (needed to conserve tool name and version)
  190 + this->GApplication::clear();
  191 +
  192 + // Initialise members
  193 + this->ctool::init_members();
  194 + this->ctobservation::init_members();
  195 + init_members();
  196 +
  197 + // Write header into logger
  198 + log_header();
  199 +
  200 + // Return
  201 + return;
  202 +}
  203 +
  204 +
  205 +/***********************************************************************//**
  206 + * @brief Select event data
  207 + *
  208 + * This method reads in the application parameters and loops over all
  209 + * observations that were found to perform an event selection. Event
  210 + * selection is done by writing each observation to a temporary file and
  211 + * re-opening the temporary file using the cfitsio event filter syntax.
  212 + * The temporary file is deleted after this action so that no disk overflow
  213 + * will occur.
  214 + ***************************************************************************/
  215 +void ctphase::run(void)
  216 +{
  217 + // Switch screen logging on in debug mode
  218 + if (logDebug()) {
  219 + log.cout(true);
  220 + }
  221 +
  222 + // Get parameters
  223 + get_parameters();
  224 +
  225 + // Write input observation container into logger
  226 + log_observations(NORMAL, m_obs, "Input observation");
  227 +
  228 + // Write header into logger
  229 + log_header1(TERSE, "Event selection");
  230 +
  231 + // Initialise counters
  232 + int n_observations = 0;
  233 +
  234 + // Loop over all observation in the container
  235 + for (int i = 0; i < m_obs.size(); ++i) {
  236 +
  237 + // Write header for the current observation
  238 + log_header3(TERSE, get_obs_header(m_obs[i]));
  239 +
  240 + // Initialise event input and output filenames and the event
  241 + // and GTI extension names
  242 + m_infiles.push_back("");
  243 + m_evtname.push_back(gammalib::extname_cta_events);
  244 + m_gtiname.push_back(gammalib::extname_gti);
  245 +
  246 + // Get CTA observation
  247 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  248 +
  249 + // Skip observation if it's not CTA
  250 + if (obs == NULL) {
  251 + std::string msg = " Skipping "+m_obs[i]->instrument()+
  252 + " observation";
  253 + log_string(NORMAL, msg);
  254 + continue;
  255 + }
  256 +
  257 + // Skip observation if we have a binned observation
  258 + if (obs->eventtype() == "CountsCube") {
  259 + std::string msg = " Skipping binned "+m_obs[i]->instrument()+
  260 + " observation";
  261 + log_string(NORMAL, msg);
  262 + continue;
  263 + }
  264 +
  265 + // Increment counter
  266 + n_observations++;
  267 +
  268 + // Save event file name (for possible saving)
  269 + m_infiles[i] = obs->eventfile();
  270 +
  271 + // Extract event and GTI extension names from input FITS file
  272 + GFilename fname(m_infiles[i]);
  273 + if (fname.has_extname()) {
  274 + m_evtname[i] = fname.extname();
  275 + }
  276 + m_gtiname[i] = get_gtiname(fname.url(), m_evtname[i]);
  277 +
  278 + // Write input file information into logger
  279 + log_value(NORMAL, "Input filename", m_infiles[i]);
  280 + log_value(NORMAL, "Event extension name", m_evtname[i]);
  281 + log_value(NORMAL, "GTI extension name", m_gtiname[i]);
  282 +
  283 + // Fall through in case that the event file is empty
  284 + if (obs->events()->size() == 0) {
  285 + log_string(NORMAL, " Warning: No events in event file \""+
  286 + m_infiles[i]+"\". Event selection skipped.");
  287 + continue;
  288 + }
  289 +
  290 + // If we have an input file then check it
  291 + if (!m_infiles[i].empty()) {
  292 + std::string message = check_infile(m_infiles[i], m_evtname[i]);
  293 + if (!message.empty()) {
  294 + throw GException::invalid_value(G_RUN, message);
  295 + }
  296 + }
  297 +
  298 + // Get temporary file name
  299 + #if G_USE_MKSTEMP
  300 + char tpl[] = "ctphaseXXXXXX";
  301 + int fileid = mkstemp(tpl);
  302 + std::string filename(tpl);
  303 + #else
  304 + std::string filename = std::tmpnam(NULL);
  305 + #endif
  306 +
  307 + // Save observation in temporary file. We add here the events and
  308 + // GTI extension name so that the GCTAObservation::save method can
  309 + // use this information for writing the proper extension names into
  310 + // the temporary file
  311 + obs->save(filename+"["+m_evtname[i]+";"+ m_gtiname[i]+"]", true);
  312 +
  313 + // Log saved FITS file.
  314 + if (logExplicit()) {
  315 + GFits tmpfile(filename);
  316 + log.header3("FITS file content of temporary file");
  317 + log << tmpfile << std::endl;
  318 + tmpfile.close();
  319 + }
  320 +
  321 + // If we have a temporary file then check it
  322 + if (!filename.empty()) {
  323 + std::string message = check_infile(filename, m_evtname[i]);
  324 + if (!message.empty()) {
  325 + throw GException::invalid_value(G_RUN, message);
  326 + }
  327 + }
  328 +
  329 + // Load observation from temporary file, including event selection
  330 + phase_events(obs, filename, m_evtname[i], m_gtiname[i]);
  331 +
  332 + // Close temporary file
  333 + #if G_USE_MKSTEMP
  334 + close(fileid);
  335 + #endif
  336 +
  337 + // Remove temporary file
  338 + std::remove(filename.c_str());
  339 +
  340 + } // endfor: looped over all observations
  341 +
  342 + // If more than a single observation has been handled then make sure that
  343 + // an XML file will be used for storage
  344 + if (n_observations > 1) {
  345 + m_use_xml = true;
  346 + }
  347 +
  348 + // Write observation(s) into logger
  349 + log_observations(NORMAL, m_obs, "Output observation");
  350 +
  351 + // Optionally publish event list(s)
  352 + if ((*this)["publish"].boolean()) {
  353 + publish();
  354 + }
  355 +
  356 + // Return
  357 + return;
  358 +}
  359 +
  360 +
  361 +/***********************************************************************//**
  362 + * @brief Save the selected event list(s)
  363 + *
  364 + * This method saves the selected event list(s) into FITS file(s). There are
  365 + * two modes, depending on the m_use_xml flag.
  366 + *
  367 + * If m_use_xml is true, all selected event list(s) will be saved into FITS
  368 + * files, where the output filenames are constructued from the input
  369 + * filenames by prepending the m_prefix string to name. Any path information
  370 + * will be stripped form the input name, hence event files will be written
  371 + * into the local working directory (unless some path information is present
  372 + * in the prefix). In addition, an XML file will be created that gathers
  373 + * the filename information for the selected event list(s). If an XML file
  374 + * was present on input, all metadata information will be copied from this
  375 + * input file.
  376 + *
  377 + * If m_use_xml is false, the selected event list will be saved into a FITS
  378 + * file.
  379 + ***************************************************************************/
  380 +void ctphase::save(void)
  381 +{
  382 + // Write header into logger
  383 + log_header1(TERSE, gammalib::number("Save event list", m_obs.size()));
  384 +
  385 + // Case A: Save event file(s) and XML metadata information
  386 + if (m_use_xml) {
  387 + save_xml();
  388 + }
  389 +
  390 + // Case B: Save event file as FITS file
  391 + else {
  392 + save_fits();
  393 + }
  394 +
  395 + // Return
  396 + return;
  397 +}
  398 +
  399 +
  400 +/***********************************************************************//**
  401 + * @brief Publish event lists
  402 + *
  403 + * @param[in] name Event list name.
  404 + ***************************************************************************/
  405 +void ctphase::publish(const std::string& name)
  406 +{
  407 + // Write header into logger
  408 + log_header1(TERSE, gammalib::number("Publish event list", m_obs.size()));
  409 +
  410 + // Loop over all observation in the container
  411 + for (int i = 0; i < m_obs.size(); ++i) {
  412 +
  413 + // Get CTA observation
  414 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  415 +
  416 + // Handle only CTA observations
  417 + if (obs != NULL) {
  418 +
  419 + // Continue only if there is an event list
  420 + if (obs->events()->size() != 0) {
  421 +
  422 + // Set default name if user name is empty
  423 + std::string user_name(name);
  424 + if (user_name.empty()) {
  425 + user_name = CTPHASE_NAME;
  426 + }
  427 +
  428 + // If there are several event lists then add an index
  429 + if (m_use_xml) {
  430 + user_name += gammalib::str(i);
  431 + }
  432 +
  433 + // Write event list name into logger
  434 + log_value(NORMAL, "Event list name", user_name);
  435 +
  436 + // Write events into in-memory FITS file
  437 + GFits fits;
  438 + obs->write(fits);
  439 +
  440 + // Publish
  441 + fits.publish(gammalib::extname_cta_events, user_name);
  442 +
  443 + } // endif: there were events
  444 +
  445 + } // endif: observation was a CTA observation
  446 +
  447 + } // endfor: looped over observations
  448 +
  449 + // Return
  450 + return;
  451 +}
  452 +
  453 +
  454 +/*==========================================================================
  455 + = =
  456 + = Private methods =
  457 + = =
  458 + ==========================================================================*/
  459 +
  460 +/***********************************************************************//**
  461 + * @brief Initialise class members
  462 + ***************************************************************************/
  463 +void ctphase::init_members(void)
  464 +{
  465 + // Initialise parameters
  466 + m_outobs.clear();
  467 + m_prefix.clear();
  468 + m_phase.clear();
  469 + m_chatter = static_cast<GChatter>(2);
  470 +
  471 + // Initialise protected members
  472 + m_infiles.clear();
  473 + m_evtname.clear();
  474 + m_gtiname.clear();
  475 +
  476 + // Return
  477 + return;
  478 +}
  479 +
  480 +
  481 +/***********************************************************************//**
  482 + * @brief Copy class members
  483 + *
  484 + * @param[in] app Application.
  485 + ***************************************************************************/
  486 +void ctphase::copy_members(const ctphase& app)
  487 +{
  488 + // Copy parameters
  489 + m_outobs = app.m_outobs;
  490 + m_prefix = app.m_prefix;
  491 + m_chatter = app.m_chatter;
  492 +
  493 + // Copy protected members
  494 + m_infiles = app.m_infiles;
  495 + m_evtname = app.m_evtname;
  496 + m_gtiname = app.m_gtiname;
  497 + m_phase = app.m_phase;
  498 +
  499 + // Return
  500 + return;
  501 +}
  502 +
  503 +
  504 +/***********************************************************************//**
  505 + * @brief Delete class members
  506 + ***************************************************************************/
  507 +void ctphase::free_members(void)
  508 +{
  509 + // Return
  510 + return;
  511 +}
  512 +
  513 +
  514 +/***********************************************************************//**
  515 + * @brief Get application parameters
  516 + *
  517 + * Get all user parameters from parameter file or (if required) by querying
  518 + * the user. Times are assumed to be in the native CTA MJD format.
  519 + *
  520 + * This method also loads observations if no observations are yet allocated.
  521 + * Observations are either loaded from a single CTA even list, or from a
  522 + * XML file using the metadata information that is stored in that file.
  523 + ***************************************************************************/
  524 +void ctphase::get_parameters(void)
  525 +{
  526 +
  527 + // Setup observations from "inobs" parameter. Do not request response
  528 + // information and do not accept counts cubes.
  529 + setup_observations(m_obs, false, true, false);
  530 +
  531 + // Check for sanity of energy selection parameters
  532 + if ((*this)["p0"].is_valid() &&
  533 + (*this)["mjd"].is_valid() &&
  534 + (*this)["f0"].is_valid()) {
  535 +
  536 + m_phase = GModelTemporalPhaseCurve();
  537 + m_phase.norm(1.0);
  538 + m_phase.mjd((*this)["mjd"].real());
  539 + m_phase.phase((*this)["p0"].real());
  540 + m_phase.f0((*this)["f0"].real());
  541 + m_phase.f1((*this)["f1"].real());
  542 + m_phase.f2((*this)["f2"].real());
  543 +
  544 + }
  545 + else {
  546 + GException::invalid_value("ctphase::get_parameters()",
  547 + "Invalid parameter value for p0, mjd, or f0.");
  548 + }
  549 +
  550 + // Get other User parameters
  551 + m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
  552 +
  553 + // Optionally read ahead parameters so that they get correctly
  554 + // dumped into the log file
  555 + if (read_ahead()) {
  556 + m_outobs = (*this)["outobs"].filename();
  557 + m_prefix = (*this)["prefix"].string();
  558 + }
  559 +
  560 + // Write parameters into logger
  561 + log_parameters(TERSE);
  562 +
  563 + // Return
  564 + return;
  565 +}
  566 +
  567 +
  568 +/***********************************************************************//**
  569 + * @brief Define phase for events
  570 + *
  571 + * @param[in,out] obs CTA observation.
  572 + * @param[in] filename File name.
  573 + * @param[in] evtname Event extension name.
  574 + * @param[in] gtiname GTI extension name.
  575 + *
  576 + * @exception GException::invalid_value
  577 + * No events extension found in FITS file.
  578 + *
  579 + * Select events from a FITS file by making use of the selection possibility
  580 + * of the cfitsio library on loading a file. A selection string is created
  581 + * from the specified criteria that is appended to the filename so that
  582 + * cfitsio will automatically filter the event data. This selection string
  583 + * is then applied when opening the FITS file. The event list in the current
  584 + * observation is replaced by selected event list read from the FITS file.
  585 + ***************************************************************************/
  586 +void ctphase::phase_events(GCTAObservation* obs,
  587 + const std::string& filename,
  588 + const std::string& evtname,
  589 + const std::string& gtiname)
  590 +{
  591 + // Write header into logger
  592 + log_header3(NORMAL, "Events phasing");
  593 +
  594 + // Initialise selection and addition strings
  595 + std::string selection;
  596 + std::string add;
  597 +
  598 + // Initialise selection flags
  599 + bool remove_all = false;
  600 +
  601 + // Get CTA event list pointer
  602 + GCTAEventList* list =
  603 + static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
  604 +
  605 + for (int evnt=0; evnt<list->size(); evnt++) {
  606 + // Get the next event
  607 + GCTAEventAtom* event = (*list)[evnt];
  608 +
  609 + // Get the event time
  610 + GTime time = event->time();
  611 +
  612 + // Set the event phase
  613 + event->phase(m_phase.phase(time));
  614 + }
  615 +
  616 + // Make sure the file knows to save the phase information
  617 + if (!list->has_phase()) {
  618 + list->has_phase(true);
  619 + }
  620 +/*
  621 + // Save GTI for later usage
  622 + GGti gti = list->gti();
  623 +
  624 + // Build input filename including selection expression
  625 + std::string expression = filename;
  626 +
  627 + // Dump FITS filename including selection expression
  628 + log_value(NORMAL, "FITS filename", expression);
  629 +
  630 + // Open FITS file
  631 + GFits file(expression);
  632 +
  633 + // Log selected FITS file
  634 + log_header3(EXPLICIT, "FITS file content after selection");
  635 + log_string(EXPLICIT, file.print(m_chatter));
  636 +
  637 + // Check if we have an events HDU
  638 + if (!file.contains(evtname)) {
  639 + std::string msg = "No events extension \""+evtname+"\" found in "
  640 + "FITS file. The expression \""+expression+"\" "
  641 + "was used to open the FITS file.";
  642 + throw GException::invalid_value(G_SELECT_EVENTS, msg);
  643 + }
  644 +
  645 + // Determine number of events in the events HDU. If removal of all events
  646 + // has been requested then set the number of events to zero.
  647 + int nevents = (remove_all) ? 0 : file.table(evtname)->nrows();
  648 +
  649 + // If the selected event list is empty then append an empty event list
  650 + // to the observation
  651 + if (nevents < 1) {
  652 +
  653 + // Create empty event list
  654 + GCTAEventList eventlist;
  655 +
  656 + // Append list to observation
  657 + obs->events(eventlist);
  658 +
  659 + }
  660 +
  661 + // ... otherwise load the data from the temporary file
  662 + else {
  663 + obs->read(file);
  664 + }
  665 +
  666 + // Get CTA event list pointer
  667 + list = static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
  668 +
  669 + // Make sure that events are fetched since the temporary file will be
  670 + // closed later
  671 + list->fetch();
  672 +
  673 + // Set event list GTI (in any case as any event list has a GTI)
  674 + list->gti(gti);
  675 +*/
  676 + // Return
  677 + return;
  678 +}
  679 +
  680 +
  681 +/***********************************************************************//**
  682 + * @brief Check input filename
  683 + *
  684 + * @param[in] filename File name.
  685 + * @param[in] evtname Event extension name.
  686 + *
  687 + * This method checks if the input FITS file is correct.
  688 + ***************************************************************************/
  689 +std::string ctphase::check_infile(const std::string& filename,
  690 + const std::string& evtname) const
  691 +{
  692 + // Initialise message string
  693 + std::string message = "";
  694 +
  695 + // Open FITS file
  696 + GFits fits(filename);
  697 +
  698 + // Check for existence of events extensions
  699 + if (!fits.contains(evtname)) {
  700 + message = "No \""+evtname+"\" extension found in input file \""+
  701 + filename + "\".";
  702 + }
  703 +
  704 + // ... otherwise check column names
  705 + else {
  706 +
  707 + // Get pointer to FITS table
  708 + GFitsTable* table = fits.table(evtname);
  709 +
  710 + // Initialise list of missing columns
  711 + std::vector<std::string> missing;
  712 +
  713 + // Check for existence of TIME column
  714 + if (!table->contains("TIME")) {
  715 + missing.push_back("TIME");
  716 + }
  717 +
  718 + // Check for existence of ENERGY column
  719 + if (!table->contains("ENERGY")) {
  720 + missing.push_back("ENERGY");
  721 + }
  722 +
  723 + // Check for existence of RA column
  724 + if (!table->contains("RA")) {
  725 + missing.push_back("RA");
  726 + }
  727 +
  728 + // Check for existence of DEC column
  729 + if (!table->contains("DEC")) {
  730 + missing.push_back("DEC");
  731 + }
  732 +
  733 + // Set error message for missing columns
  734 + if (!missing.empty()) {
  735 + message = "The following columns are missing in the "
  736 + "\""+evtname+"\" extension of input file \""+
  737 + filename + "\": ";
  738 + for (int i = 0; i < missing.size(); ++i) {
  739 + message += "\"" + missing[i] + "\"";
  740 + if (i < missing.size()-1) {
  741 + message += ", ";
  742 + }
  743 + }
  744 + }
  745 +
  746 + } // endelse: checked column names
  747 +
  748 + // Return
  749 + return message;
  750 +}
  751 +
  752 +
  753 +/***********************************************************************//**
  754 + * @brief Set output file name.
  755 + *
  756 + * @param[in] filename Input file name.
  757 + *
  758 + * Converts an input file name into an output filename by prepending the
  759 + * prefix stored in the member m_prefix to the input file name. Any path as
  760 + * well as extension will be stripped from the input file name. Also a
  761 + * trailing ".gz" will be stripped as one cannot write into gzipped files.
  762 + ***************************************************************************/
  763 +std::string ctphase::set_outfile_name(const std::string& filename) const
  764 +{
  765 + // Create filename
  766 + GFilename fname(filename);
  767 +
  768 + // Split input filename without any extensions into path elements
  769 + std::vector<std::string> elements = gammalib::split(fname.url(), "/");
  770 +
  771 + // The last path element is the filename
  772 + std::string outname = m_prefix + elements[elements.size()-1];
  773 +
  774 + // Strip any ".gz"
  775 + outname = gammalib::strip_chars(outname, ".gz");
  776 +
  777 + // Return output filename
  778 + return outname;
  779 +}
  780 +
  781 +
  782 +/***********************************************************************//**
  783 + * @brief Get Good Time Intervals extension name
  784 + *
  785 + * @param[in] filename Input file name.
  786 + * @param[in] evtname Events extension name.
  787 + *
  788 + * Extracts the Good Time Intervals extension name from the event file. We
  789 + * do this by loading the events and accessing the Good Time Intervals
  790 + * extension name using the GCTAEventList::gtiname() method. If the file name
  791 + * is empty, the method returns `GTI`.
  792 + ***************************************************************************/
  793 +std::string ctphase::get_gtiname(const std::string& filename,
  794 + const std::string& evtname) const
  795 +{
  796 + // Initialise GTI name
  797 + std::string gtiname = gammalib::extname_gti;
  798 +
  799 + // Continue only if the filename is not empty
  800 + if (!filename.empty()) {
  801 +
  802 + // Load events
  803 + GCTAEventList events(filename+"["+evtname+"]");
  804 +
  805 + // Get GTI name
  806 + gtiname = events.gtiname();
  807 +
  808 + }
  809 +
  810 + // Return GTI name
  811 + return (gtiname);
  812 +}
  813 +
  814 +
  815 +/***********************************************************************//**
  816 + * @brief Save event list in FITS format.
  817 + *
  818 + * Save the event list as a FITS file. The file name of the FITS file is
  819 + * specified by the "outobs" parameter.
  820 + ***************************************************************************/
  821 +void ctphase::save_fits(void)
  822 +{
  823 + // Save only if there are observations
  824 + if (m_obs.size() > 0) {
  825 +
  826 + // Get output filename
  827 + m_outobs = (*this)["outobs"].filename();
  828 +
  829 + // Get CTA observation from observation container
  830 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
  831 +
  832 + // Save only if it's a CTA observation
  833 + if (obs != NULL) {
  834 +
  835 + // Save only if file name is non-empty
  836 + if (m_infiles[0].length() > 0) {
  837 +
  838 + // Create file name object
  839 + GFilename fname(m_outobs);
  840 +
  841 + // Extract filename and event extension name
  842 + std::string outfile = fname.url();
  843 +
  844 + // Append event extension name. We handle here the possibility
  845 + // to write the events into a different extension.
  846 + if (fname.has_extname()) {
  847 + outfile += "["+fname.extname()+"]";
  848 + }
  849 + else {
  850 + outfile += "["+m_evtname[0]+"]";
  851 + }
  852 +
  853 + // Log filename
  854 + log_value(NORMAL, "Event list file", outfile);
  855 +
  856 + // Save event list
  857 + save_event_list(obs, m_infiles[0], m_evtname[0], m_gtiname[0],
  858 + outfile);
  859 +
  860 + } // endif: filename was non empty
  861 +
  862 + } // endif: observation was CTA observation
  863 +
  864 + } // endif: there were observations
  865 +
  866 + // Return
  867 + return;
  868 +}
  869 +
  870 +
  871 +/***********************************************************************//**
  872 + * @brief Save event list(s) in XML format.
  873 + *
  874 + * Save the event list(s) into FITS files and write the file path information
  875 + * into a XML file. The filename of the XML file is specified by the outfile
  876 + * parameter, the filename(s) of the event lists are built by prepending a
  877 + * prefix to the input event list filenames. Any path present in the input
  878 + * filename will be stripped, i.e. the event list(s) will be written in the
  879 + * local working directory (unless a path is specified in the prefix).
  880 + ***************************************************************************/
  881 +void ctphase::save_xml(void)
  882 +{
  883 + // Get output filename and prefix
  884 + m_outobs = (*this)["outobs"].filename();
  885 + m_prefix = (*this)["prefix"].string();
  886 +
  887 + // Issue warning if output filename has no .xml suffix
  888 + log_string(TERSE, warn_xml_suffix(m_outobs));
  889 +
  890 + // Loop over all observation in the container
  891 + for (int i = 0; i < m_obs.size(); ++i) {
  892 +
  893 + // Get CTA observation
  894 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  895 +
  896 + // Skip observations that are no CTA observations
  897 + if (obs == NULL) {
  898 + continue;
  899 + }
  900 +
  901 + // Skip observations that have empty names
  902 + if (m_infiles[i].length() == 0) {
  903 + continue;
  904 + }
  905 +
  906 + // Set event output file name
  907 + std::string outfile = set_outfile_name(m_infiles[i]);
  908 +
  909 + // Append event extension name
  910 + outfile += "["+m_evtname[i]+"]";
  911 +
  912 + // Log filename
  913 + log_value(NORMAL, "Event list file", outfile);
  914 +
  915 + // Store output file name in observation
  916 + obs->eventfile(outfile);
  917 +
  918 + // Save event list
  919 + save_event_list(obs, m_infiles[i], m_evtname[i], m_gtiname[i],
  920 + outfile);
  921 +
  922 + } // endfor: looped over observations
  923 +
  924 + // Write observation definition XML file name into logger
  925 + log_value(NORMAL, "Obs. definition file", m_outobs);
  926 +
  927 + // Save observations in XML file
  928 + m_obs.save(m_outobs);
  929 +
  930 + // Return
  931 + return;
  932 +}
  933 +
  934 +
  935 +/***********************************************************************//**
  936 + * @brief Save event list into FITS file
  937 + *
  938 + * @param[in] obs Pointer to CTA observation.
  939 + * @param[in] infile Input file name.
  940 + * @param[in] evtname Event extension name.
  941 + * @param[in] gtiname GTI extension name.
  942 + * @param[in] outfile Output file name.
  943 + *
  944 + * Saves an event list and the corresponding Good Time Intervals into a FITS
  945 + * file and copy all others extensions from the input file to the output
  946 + * file.
  947 + *
  948 + * If an extension name is specified in the @p outfile argument, the events
  949 + * and eventually also the Good Time Intervals will be extracted from the
  950 + * argument and used for writing the events. The format is
  951 + *
  952 + * <filename>[<event extension name;GTI extension name>]
  953 + *
  954 + * where <filename> needs to be replaced by the name of the FITS file,
  955 + * and <event extension name;GTI extension name> by the name of the events
  956 + * and Good Time Intervals extensions. For example
  957 + *
  958 + * myfits.fits[EVENTS1;GTI1]
  959 + *
  960 + * will write the selected events into the "EVENTS1" extension and the
  961 + * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS
  962 + * file. If the Good Time Intervals extension name is skipped, e.g.
  963 + *
  964 + * myfits.fits[EVENTS1]
  965 + *
  966 + * the original extension name for the Good Time Intervals will be kept.
  967 + * Analogously, only the Good Time Intervals extension name can be changed
  968 + * by specifying
  969 + *
  970 + * myfits.fits[;GTI1]
  971 + *
  972 + * In none of the cases will the original events and Good Time Intervals be
  973 + * copied over to the output file.
  974 + ***************************************************************************/
  975 +void ctphase::save_event_list(const GCTAObservation* obs,
  976 + const std::string& infile,
  977 + const std::string& evtname,
  978 + const std::string& gtiname,
  979 + const std::string& outfile) const
  980 +{
  981 + // Save only if we have an event list
  982 + if (obs->eventtype() == "EventList") {
  983 +
  984 + // Set output FITS file event extension names
  985 + GFilename outname(outfile);
  986 + std::string outevt = evtname;
  987 + std::string outgti = gtiname;
  988 + if (outname.has_extname()) {
  989 + std::vector<std::string> extnames =
  990 + gammalib::split(outname.extname(), ";");
  991 + if (extnames.size() > 0) {
  992 + std::string extname = gammalib::strip_whitespace(extnames[0]);
  993 + if (!extname.empty()) {
  994 + outevt = extname;
  995 + }
  996 + }
  997 + if (extnames.size() > 1) {
  998 + std::string extname = gammalib::strip_whitespace(extnames[1]);
  999 + if (!extname.empty()) {
  1000 + outgti = extname;
  1001 + }
  1002 + }
  1003 + }
  1004 +
  1005 + // Create output FITS file
  1006 + GFits outfits;
  1007 +
  1008 + // Write observation into FITS file
  1009 + obs->write(outfits, outevt, outgti);
  1010 +
  1011 + // Copy all extensions other than evtname and gtiname extensions
  1012 + // from the input to the output event list. The evtname and
  1013 + // gtiname extensions are written by the save method, all others
  1014 + // that may eventually be present have to be copied over
  1015 + // explicitly.
  1016 + GFits infits(infile);
  1017 + for (int extno = 1; extno < infits.size(); ++extno) {
  1018 + GFitsHDU* hdu = infits.at(extno);
  1019 + if (hdu->extname() != evtname &&
  1020 + hdu->extname() != gtiname &&
  1021 + hdu->extname() != outevt &&
  1022 + hdu->extname() != outgti) {
  1023 + outfits.append(*hdu);
  1024 + }
  1025 + }
  1026 +
  1027 + // Close input file
  1028 + infits.close();
  1029 +
  1030 + // Save file to disk and close it (we need both operations)
  1031 + outfits.saveto(outname.url(), clobber());
  1032 + outfits.close();
  1033 +
  1034 + } // endif: observation was unbinned
  1035 +
  1036 + // Return
  1037 + return;
  1038 +}
... ...
src/ctphase/ctphase.hpp 0 → 100644
  View file @6f2f518
  1 +/***************************************************************************
  2 + * ctphase - Data selection tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 by Leonardo Di Venere *
  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 ctphase.hpp
  23 + * @brief Data selection tool definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef CTPHASE_HPP
  28 +#define CTPHASE_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <vector>
  32 +#include <string>
  33 +#include "GammaLib.hpp"
  34 +#include "GCTALib.hpp"
  35 +#include "ctobservation.hpp"
  36 +
  37 +/* __Definitions _________________________________________________________ */
  38 +#define CTPHASE_NAME "ctphase"
  39 +#define CTPHASE_VERSION "0.0.1"
  40 +
  41 +
  42 +/***********************************************************************//**
  43 + * @class ctphase
  44 + *
  45 + * @brief Append phase columm to observation files
  46 + ***************************************************************************/
  47 +class ctphase : public ctobservation {
  48 +
  49 +public:
  50 + // Constructors and destructors
  51 + ctphase(void);
  52 + explicit ctphase(const GObservations& obs);
  53 + ctphase(int argc, char *argv[]);
  54 + ctphase(const ctphase& app);
  55 + virtual ~ctphase(void);
  56 +
  57 + // Operators
  58 + ctphase& operator=(const ctphase& app);
  59 +
  60 + // Methods
  61 + void clear(void);
  62 + void run(void);
  63 + void save(void);
  64 + void publish(const std::string& name = "");
  65 +
  66 + // TODO: Add Methods for updating model parameters
  67 +
  68 +protected:
  69 + // Protected methods
  70 + void init_members(void);
  71 + void copy_members(const ctphase& app);
  72 + void free_members(void);
  73 + void get_parameters(void);
  74 + void phase_events(GCTAObservation* obs,
  75 + const std::string& filename,
  76 + const std::string& evtname,
  77 + const std::string& gtiname);
  78 + std::string check_infile(const std::string& filename,
  79 + const std::string& evtname) const;
  80 + std::string set_outfile_name(const std::string& filename) const;
  81 + std::string get_gtiname(const std::string& filename,
  82 + const std::string& evtname) const;
  83 + void save_fits(void);
  84 + void save_xml(void);
  85 + void save_event_list(const GCTAObservation* obs,
  86 + const std::string& infile,
  87 + const std::string& evtname,
  88 + const std::string& gtiname,
  89 + const std::string& outfile) const;
  90 +
  91 + // User parameters
  92 + std::string m_outobs; //!< Output event list or XML file
  93 + std::string m_prefix; //!< Prefix for multiple event lists
  94 + GModelTemporalPhaseCurve m_phase; //!< Phase information object
  95 + GChatter m_chatter; //!< Chattiness
  96 +
  97 + // Protected members
  98 + std::vector<std::string> m_infiles; //!< Input event filenames
  99 + std::vector<std::string> m_evtname; //!< Event extension names
  100 + std::vector<std::string> m_gtiname; //!< GTI extension names
  101 +};
  102 +
  103 +
  104 +#endif /* CTPHASE_HPP */
... ...
src/ctphase/ctphase.par 0 → 100644
  View file @6f2f518
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: ctphase.par #
  8 +# Version: 0.0.1 #
  9 +# #
  10 +# Author: Leonardo Di Venere #
  11 +# leonardo.divenere@ba.infn.it #
  12 +# INFN #
  13 +# #
  14 +# Purpose: Parameter file for the ctphase application #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# General parameters
  19 +#===================
  20 +inobs, f, a, "events.fits",,, "Input event list or observation definition XML file"
  21 +outobs, f, a, "phased_events.fits",,, "Output event list or observation definition XML file"
  22 +prefix, s, h, "phased_",,, "Prefix for event lists in observation definition XML file"
  23 +
  24 +#
  25 +# Selection parameters
  26 +#=====================
  27 +mjd, r, a, 52400,,, "Reference MJD for phase information"
  28 +f0, r, a, 1,,, "Frequency at reference MJD"
  29 +f1, r, a, 0,,, "First freq. derivative at reference MJD"
  30 +f2, r, a, 0,,, "Second freq. derivative at reference MJD"
  31 +p0, r, a, 0.1,0.0,1.0, "Phase at reference MJD"
  32 +
  33 +#
  34 +# Standard parameters
  35 +#====================
  36 +publish, b, h, no,,, "Publish event list on VO Hub?"
  37 +chatter, i, h, 2,0,4, "Chattiness of output"
  38 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  39 +debug, b, h, no,,, "Debugging mode activated"
  40 +mode, s, h, "h",,, "Mode of automatic parameters"
  41 +logfile, f, h, "ctphase.log",,, "Log filename"
... ...
src/ctphase/main.cpp 0 → 100644
  View file @6f2f518
  1 +/***************************************************************************
  2 + * ctphase - Append phase information to CTA events file *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 by Leonardo Di Venere *
  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 ctphase/main.cpp
  23 + * @brief Append phase information to CTA events file
  24 + * @author Leonardo Di Venere
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "support.hpp"
  32 +#include "ctphase.hpp"
  33 +
  34 +
  35 +/***********************************************************************//**
  36 + * @brief Main entry point
  37 + *
  38 + * @param[in] argc Number of arguments
  39 + * @param[in] argv Arguments
  40 + *
  41 + * This is the main entry point of the ctphase application. It allocates a
  42 + * ctphase object and executes the application. Any exceptions that occur
  43 + * will be catched and corresponding error messages written in the
  44 + * application logger and into the standard output.
  45 + ***************************************************************************/
  46 +int main (int argc, char *argv[])
  47 +{
  48 + // Initialise return code
  49 + int rc = 1;
  50 +
  51 + // Initialise pointer on application
  52 + ctphase* application = NULL;
  53 +
  54 + // Execute application
  55 + try {
  56 +
  57 + // Create instance of application
  58 + application = new ctphase(argc, argv);
  59 +
  60 + // Execute ctool
  61 + rc = execute_ctool(application);
  62 +
  63 + }
  64 +
  65 + catch (std::exception &e) {
  66 +
  67 + // Report exception
  68 + report_ctool_failure("ctphase", e.what());
  69 +
  70 + }
  71 +
  72 + // Delete application
  73 + delete application;
  74 +
  75 + // Return
  76 + return rc;
  77 +}
... ...