From 6f2f518db3dac89ae204497c99a0c9c7e147a828 Mon Sep 17 00:00:00 2001 From: Josh Cardenzana Date: Tue, 4 Apr 2017 14:26:49 +0200 Subject: [PATCH] Updated ctphase tool (#1982) --- configure.ac | 1 + src/Makefile.am | 7 +++++-- src/ctphase/Makefile.am | 47 +++++++++++++++++++++++++++++++++++++++++++++++ src/ctphase/ctphase.cpp | 1038 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctphase/ctphase.hpp | 104 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctphase/ctphase.par | 41 +++++++++++++++++++++++++++++++++++++++++ src/ctphase/main.cpp | 77 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 1313 insertions(+), 2 deletions(-) create mode 100644 src/ctphase/Makefile.am create mode 100644 src/ctphase/ctphase.cpp create mode 100644 src/ctphase/ctphase.hpp create mode 100644 src/ctphase/ctphase.par create mode 100644 src/ctphase/main.cpp diff --git a/configure.ac b/configure.ac index b16543d..99a4972 100644 --- a/configure.ac +++ b/configure.ac @@ -1018,6 +1018,7 @@ AC_CONFIG_FILES([Makefile src/ctobssim/Makefile src/ctskymap/Makefile src/ctselect/Makefile + src/ctphase/Makefile src/ctexpcube/Makefile src/ctpsfcube/Makefile src/ctedispcube/Makefile diff --git a/src/Makefile.am b/src/Makefile.am index d56177f..b30eae9 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -40,7 +40,8 @@ SUBDIRS = support \ cttsmap \ ctbutterfly \ ctulimit \ - cterror + cterror \ + ctphase # Scripts to be distributed scripts = ctools-setup \ @@ -68,7 +69,9 @@ libctools_la_LIBADD = support/libsupport.la \ cttsmap/libcttsmap.la \ ctbutterfly/libctbutterfly.la \ ctulimit/libctulimit.la \ - cterror/libcterror.la + cterror/libcterror.la \ + ctphase/libctphase.la + libctools_la_LDFLAGS = -version-info $(CTOOLS_LT_VERSION) @OPENMP_LIB_LDFLAGS@ # Distribute scripts diff --git a/src/ctphase/Makefile.am b/src/ctphase/Makefile.am new file mode 100644 index 0000000..a31e611 --- /dev/null +++ b/src/ctphase/Makefile.am @@ -0,0 +1,47 @@ +############################################################################# +# Makefile.am for ctphase tool # +# ------------------------------------------------------------------------- # +# Copyright (C) 2017 Leonardo Di Venere # +# ------------------------------------------------------------------------- # +# # +# This program is free software: you can redistribute it and/or modify # +# it under the terms of the GNU General Public License as published by # +# the Free Software Foundation, either version 3 of the License, or # +# (at your option) any later version. # +# # +# This program is distributed in the hope that it will be useful, # +# but WITHOUT ANY WARRANTY; without even the implied warranty of # +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +# GNU General Public License for more details. # +# # +# You should have received a copy of the GNU General Public License # +# along with this program. If not, see . # +# # +############################################################################# + +# Set include directories +INCLUDES = -I$(top_srcdir)/src/support + +# Program name +bin_PROGRAMS = ctphase + +# Sources +sources = ctphase.cpp ctphase.hpp + +# Libtool library (only needed for python bindings) +noinst_LTLIBRARIES = libctphase.la +libctphase_la_SOURCES = $(sources) +libctphase_la_LDFLAGS = @OPENMP_LIB_LDFLAGS@ + +# Program definition +ctphase_SOURCES = main.cpp +ctphase_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@ +ctphase_CPPFLAGS = @CPPFLAGS@ +ctphase_LDADD = libctphase.la $(top_builddir)/src/support/libsupport.la + +# Parameter files +pfilesdir = $(prefix)/syspfiles +pfiles_DATA = $(srcdir)/ctphase.par + +# Extra files to be distributed in the tar file +EXTRA_DIST = $(pfiles_DATA) diff --git a/src/ctphase/ctphase.cpp b/src/ctphase/ctphase.cpp new file mode 100644 index 0000000..e71641e --- /dev/null +++ b/src/ctphase/ctphase.cpp @@ -0,0 +1,1038 @@ +/*************************************************************************** + * ctphase - Data selection tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2010-2017 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctphase.cpp + * @brief Data selection tool definition + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctphase.hpp" +#include "GTools.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_RUN "ctphase::run()" +#define G_SELECT_EVENTS "ctphase::select_events(GCTAObservation*, "\ + "std::string&, std::string&, std::string&)" +#define G_SET_EBOUNDS "ctphase::set_ebounds(GCTAObservation*, GEbounds&)" + +/* __ Debug definitions __________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ +//#define G_USE_MKSTEMP //!< Use mkstemp for temporary filename + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +ctphase::ctphase(void) : ctobservation(CTPHASE_NAME, CTPHASE_VERSION) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Observations constructor + * + * param[in] obs Observation container. + * + * Creates an instance of the class that is initialised using the information + * provided in an observation container. + ***************************************************************************/ +ctphase::ctphase(const GObservations& obs) : + ctobservation(CTPHASE_NAME, CTPHASE_VERSION, obs) +{ + // Initialise members + init_members(); + + // Return + return; +} + + + + +/***********************************************************************//** + * @brief Command line constructor + * + * @param[in] argc Number of arguments in command line. + * @param[in] argv Array of command line arguments. + ***************************************************************************/ +ctphase::ctphase(int argc, char *argv[]) : + ctobservation(CTPHASE_NAME, CTPHASE_VERSION, argc, argv) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] app Application. + ***************************************************************************/ +ctphase::ctphase(const ctphase& app) : ctobservation(app) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +ctphase::~ctphase(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] app Application. + * @return Application. + ***************************************************************************/ +ctphase& ctphase::operator=(const ctphase& app) +{ + // Execute only if object is not identical + if (this != &app) { + + // Copy base class members + this->ctobservation::operator=(app); + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + } // endif: object was not identical + + // Return this object + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear ctphase tool + * + * Clears ctphase tool. + ***************************************************************************/ +void ctphase::clear(void) +{ + // Free members + free_members(); + this->ctobservation::free_members(); + this->ctool::free_members(); + + // Clear base class (needed to conserve tool name and version) + this->GApplication::clear(); + + // Initialise members + this->ctool::init_members(); + this->ctobservation::init_members(); + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Select event data + * + * This method reads in the application parameters and loops over all + * observations that were found to perform an event selection. Event + * selection is done by writing each observation to a temporary file and + * re-opening the temporary file using the cfitsio event filter syntax. + * The temporary file is deleted after this action so that no disk overflow + * will occur. + ***************************************************************************/ +void ctphase::run(void) +{ + // Switch screen logging on in debug mode + if (logDebug()) { + log.cout(true); + } + + // Get parameters + get_parameters(); + + // Write input observation container into logger + log_observations(NORMAL, m_obs, "Input observation"); + + // Write header into logger + log_header1(TERSE, "Event selection"); + + // Initialise counters + int n_observations = 0; + + // Loop over all observation in the container + for (int i = 0; i < m_obs.size(); ++i) { + + // Write header for the current observation + log_header3(TERSE, get_obs_header(m_obs[i])); + + // Initialise event input and output filenames and the event + // and GTI extension names + m_infiles.push_back(""); + m_evtname.push_back(gammalib::extname_cta_events); + m_gtiname.push_back(gammalib::extname_gti); + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[i]); + + // Skip observation if it's not CTA + if (obs == NULL) { + std::string msg = " Skipping "+m_obs[i]->instrument()+ + " observation"; + log_string(NORMAL, msg); + continue; + } + + // Skip observation if we have a binned observation + if (obs->eventtype() == "CountsCube") { + std::string msg = " Skipping binned "+m_obs[i]->instrument()+ + " observation"; + log_string(NORMAL, msg); + continue; + } + + // Increment counter + n_observations++; + + // Save event file name (for possible saving) + m_infiles[i] = obs->eventfile(); + + // Extract event and GTI extension names from input FITS file + GFilename fname(m_infiles[i]); + if (fname.has_extname()) { + m_evtname[i] = fname.extname(); + } + m_gtiname[i] = get_gtiname(fname.url(), m_evtname[i]); + + // Write input file information into logger + log_value(NORMAL, "Input filename", m_infiles[i]); + log_value(NORMAL, "Event extension name", m_evtname[i]); + log_value(NORMAL, "GTI extension name", m_gtiname[i]); + + // Fall through in case that the event file is empty + if (obs->events()->size() == 0) { + log_string(NORMAL, " Warning: No events in event file \""+ + m_infiles[i]+"\". Event selection skipped."); + continue; + } + + // If we have an input file then check it + if (!m_infiles[i].empty()) { + std::string message = check_infile(m_infiles[i], m_evtname[i]); + if (!message.empty()) { + throw GException::invalid_value(G_RUN, message); + } + } + + // Get temporary file name + #if G_USE_MKSTEMP + char tpl[] = "ctphaseXXXXXX"; + int fileid = mkstemp(tpl); + std::string filename(tpl); + #else + std::string filename = std::tmpnam(NULL); + #endif + + // Save observation in temporary file. We add here the events and + // GTI extension name so that the GCTAObservation::save method can + // use this information for writing the proper extension names into + // the temporary file + obs->save(filename+"["+m_evtname[i]+";"+ m_gtiname[i]+"]", true); + + // Log saved FITS file. + if (logExplicit()) { + GFits tmpfile(filename); + log.header3("FITS file content of temporary file"); + log << tmpfile << std::endl; + tmpfile.close(); + } + + // If we have a temporary file then check it + if (!filename.empty()) { + std::string message = check_infile(filename, m_evtname[i]); + if (!message.empty()) { + throw GException::invalid_value(G_RUN, message); + } + } + + // Load observation from temporary file, including event selection + phase_events(obs, filename, m_evtname[i], m_gtiname[i]); + + // Close temporary file + #if G_USE_MKSTEMP + close(fileid); + #endif + + // Remove temporary file + std::remove(filename.c_str()); + + } // endfor: looped over all observations + + // If more than a single observation has been handled then make sure that + // an XML file will be used for storage + if (n_observations > 1) { + m_use_xml = true; + } + + // Write observation(s) into logger + log_observations(NORMAL, m_obs, "Output observation"); + + // Optionally publish event list(s) + if ((*this)["publish"].boolean()) { + publish(); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save the selected event list(s) + * + * This method saves the selected event list(s) into FITS file(s). There are + * two modes, depending on the m_use_xml flag. + * + * If m_use_xml is true, all selected event list(s) will be saved into FITS + * files, where the output filenames are constructued from the input + * filenames by prepending the m_prefix string to name. Any path information + * will be stripped form the input name, hence event files will be written + * into the local working directory (unless some path information is present + * in the prefix). In addition, an XML file will be created that gathers + * the filename information for the selected event list(s). If an XML file + * was present on input, all metadata information will be copied from this + * input file. + * + * If m_use_xml is false, the selected event list will be saved into a FITS + * file. + ***************************************************************************/ +void ctphase::save(void) +{ + // Write header into logger + log_header1(TERSE, gammalib::number("Save event list", m_obs.size())); + + // Case A: Save event file(s) and XML metadata information + if (m_use_xml) { + save_xml(); + } + + // Case B: Save event file as FITS file + else { + save_fits(); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Publish event lists + * + * @param[in] name Event list name. + ***************************************************************************/ +void ctphase::publish(const std::string& name) +{ + // Write header into logger + log_header1(TERSE, gammalib::number("Publish event list", m_obs.size())); + + // Loop over all observation in the container + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[i]); + + // Handle only CTA observations + if (obs != NULL) { + + // Continue only if there is an event list + if (obs->events()->size() != 0) { + + // Set default name if user name is empty + std::string user_name(name); + if (user_name.empty()) { + user_name = CTPHASE_NAME; + } + + // If there are several event lists then add an index + if (m_use_xml) { + user_name += gammalib::str(i); + } + + // Write event list name into logger + log_value(NORMAL, "Event list name", user_name); + + // Write events into in-memory FITS file + GFits fits; + obs->write(fits); + + // Publish + fits.publish(gammalib::extname_cta_events, user_name); + + } // endif: there were events + + } // endif: observation was a CTA observation + + } // endfor: looped over observations + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctphase::init_members(void) +{ + // Initialise parameters + m_outobs.clear(); + m_prefix.clear(); + m_phase.clear(); + m_chatter = static_cast(2); + + // Initialise protected members + m_infiles.clear(); + m_evtname.clear(); + m_gtiname.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctphase::copy_members(const ctphase& app) +{ + // Copy parameters + m_outobs = app.m_outobs; + m_prefix = app.m_prefix; + m_chatter = app.m_chatter; + + // Copy protected members + m_infiles = app.m_infiles; + m_evtname = app.m_evtname; + m_gtiname = app.m_gtiname; + m_phase = app.m_phase; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctphase::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Get application parameters + * + * Get all user parameters from parameter file or (if required) by querying + * the user. Times are assumed to be in the native CTA MJD format. + * + * This method also loads observations if no observations are yet allocated. + * Observations are either loaded from a single CTA even list, or from a + * XML file using the metadata information that is stored in that file. + ***************************************************************************/ +void ctphase::get_parameters(void) +{ + + // Setup observations from "inobs" parameter. Do not request response + // information and do not accept counts cubes. + setup_observations(m_obs, false, true, false); + + // Check for sanity of energy selection parameters + if ((*this)["p0"].is_valid() && + (*this)["mjd"].is_valid() && + (*this)["f0"].is_valid()) { + + m_phase = GModelTemporalPhaseCurve(); + m_phase.norm(1.0); + m_phase.mjd((*this)["mjd"].real()); + m_phase.phase((*this)["p0"].real()); + m_phase.f0((*this)["f0"].real()); + m_phase.f1((*this)["f1"].real()); + m_phase.f2((*this)["f2"].real()); + + } + else { + GException::invalid_value("ctphase::get_parameters()", + "Invalid parameter value for p0, mjd, or f0."); + } + + // Get other User parameters + m_chatter = static_cast((*this)["chatter"].integer()); + + // Optionally read ahead parameters so that they get correctly + // dumped into the log file + if (read_ahead()) { + m_outobs = (*this)["outobs"].filename(); + m_prefix = (*this)["prefix"].string(); + } + + // Write parameters into logger + log_parameters(TERSE); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Define phase for events + * + * @param[in,out] obs CTA observation. + * @param[in] filename File name. + * @param[in] evtname Event extension name. + * @param[in] gtiname GTI extension name. + * + * @exception GException::invalid_value + * No events extension found in FITS file. + * + * Select events from a FITS file by making use of the selection possibility + * of the cfitsio library on loading a file. A selection string is created + * from the specified criteria that is appended to the filename so that + * cfitsio will automatically filter the event data. This selection string + * is then applied when opening the FITS file. The event list in the current + * observation is replaced by selected event list read from the FITS file. + ***************************************************************************/ +void ctphase::phase_events(GCTAObservation* obs, + const std::string& filename, + const std::string& evtname, + const std::string& gtiname) +{ + // Write header into logger + log_header3(NORMAL, "Events phasing"); + + // Initialise selection and addition strings + std::string selection; + std::string add; + + // Initialise selection flags + bool remove_all = false; + + // Get CTA event list pointer + GCTAEventList* list = + static_cast(const_cast(obs->events())); + + for (int evnt=0; evntsize(); evnt++) { + // Get the next event + GCTAEventAtom* event = (*list)[evnt]; + + // Get the event time + GTime time = event->time(); + + // Set the event phase + event->phase(m_phase.phase(time)); + } + + // Make sure the file knows to save the phase information + if (!list->has_phase()) { + list->has_phase(true); + } +/* + // Save GTI for later usage + GGti gti = list->gti(); + + // Build input filename including selection expression + std::string expression = filename; + + // Dump FITS filename including selection expression + log_value(NORMAL, "FITS filename", expression); + + // Open FITS file + GFits file(expression); + + // Log selected FITS file + log_header3(EXPLICIT, "FITS file content after selection"); + log_string(EXPLICIT, file.print(m_chatter)); + + // Check if we have an events HDU + if (!file.contains(evtname)) { + std::string msg = "No events extension \""+evtname+"\" found in " + "FITS file. The expression \""+expression+"\" " + "was used to open the FITS file."; + throw GException::invalid_value(G_SELECT_EVENTS, msg); + } + + // Determine number of events in the events HDU. If removal of all events + // has been requested then set the number of events to zero. + int nevents = (remove_all) ? 0 : file.table(evtname)->nrows(); + + // If the selected event list is empty then append an empty event list + // to the observation + if (nevents < 1) { + + // Create empty event list + GCTAEventList eventlist; + + // Append list to observation + obs->events(eventlist); + + } + + // ... otherwise load the data from the temporary file + else { + obs->read(file); + } + + // Get CTA event list pointer + list = static_cast(const_cast(obs->events())); + + // Make sure that events are fetched since the temporary file will be + // closed later + list->fetch(); + + // Set event list GTI (in any case as any event list has a GTI) + list->gti(gti); +*/ + // Return + return; +} + + +/***********************************************************************//** + * @brief Check input filename + * + * @param[in] filename File name. + * @param[in] evtname Event extension name. + * + * This method checks if the input FITS file is correct. + ***************************************************************************/ +std::string ctphase::check_infile(const std::string& filename, + const std::string& evtname) const +{ + // Initialise message string + std::string message = ""; + + // Open FITS file + GFits fits(filename); + + // Check for existence of events extensions + if (!fits.contains(evtname)) { + message = "No \""+evtname+"\" extension found in input file \""+ + filename + "\"."; + } + + // ... otherwise check column names + else { + + // Get pointer to FITS table + GFitsTable* table = fits.table(evtname); + + // Initialise list of missing columns + std::vector missing; + + // Check for existence of TIME column + if (!table->contains("TIME")) { + missing.push_back("TIME"); + } + + // Check for existence of ENERGY column + if (!table->contains("ENERGY")) { + missing.push_back("ENERGY"); + } + + // Check for existence of RA column + if (!table->contains("RA")) { + missing.push_back("RA"); + } + + // Check for existence of DEC column + if (!table->contains("DEC")) { + missing.push_back("DEC"); + } + + // Set error message for missing columns + if (!missing.empty()) { + message = "The following columns are missing in the " + "\""+evtname+"\" extension of input file \""+ + filename + "\": "; + for (int i = 0; i < missing.size(); ++i) { + message += "\"" + missing[i] + "\""; + if (i < missing.size()-1) { + message += ", "; + } + } + } + + } // endelse: checked column names + + // Return + return message; +} + + +/***********************************************************************//** + * @brief Set output file name. + * + * @param[in] filename Input file name. + * + * Converts an input file name into an output filename by prepending the + * prefix stored in the member m_prefix to the input file name. Any path as + * well as extension will be stripped from the input file name. Also a + * trailing ".gz" will be stripped as one cannot write into gzipped files. + ***************************************************************************/ +std::string ctphase::set_outfile_name(const std::string& filename) const +{ + // Create filename + GFilename fname(filename); + + // Split input filename without any extensions into path elements + std::vector elements = gammalib::split(fname.url(), "/"); + + // The last path element is the filename + std::string outname = m_prefix + elements[elements.size()-1]; + + // Strip any ".gz" + outname = gammalib::strip_chars(outname, ".gz"); + + // Return output filename + return outname; +} + + +/***********************************************************************//** + * @brief Get Good Time Intervals extension name + * + * @param[in] filename Input file name. + * @param[in] evtname Events extension name. + * + * Extracts the Good Time Intervals extension name from the event file. We + * do this by loading the events and accessing the Good Time Intervals + * extension name using the GCTAEventList::gtiname() method. If the file name + * is empty, the method returns `GTI`. + ***************************************************************************/ +std::string ctphase::get_gtiname(const std::string& filename, + const std::string& evtname) const +{ + // Initialise GTI name + std::string gtiname = gammalib::extname_gti; + + // Continue only if the filename is not empty + if (!filename.empty()) { + + // Load events + GCTAEventList events(filename+"["+evtname+"]"); + + // Get GTI name + gtiname = events.gtiname(); + + } + + // Return GTI name + return (gtiname); +} + + +/***********************************************************************//** + * @brief Save event list in FITS format. + * + * Save the event list as a FITS file. The file name of the FITS file is + * specified by the "outobs" parameter. + ***************************************************************************/ +void ctphase::save_fits(void) +{ + // Save only if there are observations + if (m_obs.size() > 0) { + + // Get output filename + m_outobs = (*this)["outobs"].filename(); + + // Get CTA observation from observation container + GCTAObservation* obs = dynamic_cast(m_obs[0]); + + // Save only if it's a CTA observation + if (obs != NULL) { + + // Save only if file name is non-empty + if (m_infiles[0].length() > 0) { + + // Create file name object + GFilename fname(m_outobs); + + // Extract filename and event extension name + std::string outfile = fname.url(); + + // Append event extension name. We handle here the possibility + // to write the events into a different extension. + if (fname.has_extname()) { + outfile += "["+fname.extname()+"]"; + } + else { + outfile += "["+m_evtname[0]+"]"; + } + + // Log filename + log_value(NORMAL, "Event list file", outfile); + + // Save event list + save_event_list(obs, m_infiles[0], m_evtname[0], m_gtiname[0], + outfile); + + } // endif: filename was non empty + + } // endif: observation was CTA observation + + } // endif: there were observations + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save event list(s) in XML format. + * + * Save the event list(s) into FITS files and write the file path information + * into a XML file. The filename of the XML file is specified by the outfile + * parameter, the filename(s) of the event lists are built by prepending a + * prefix to the input event list filenames. Any path present in the input + * filename will be stripped, i.e. the event list(s) will be written in the + * local working directory (unless a path is specified in the prefix). + ***************************************************************************/ +void ctphase::save_xml(void) +{ + // Get output filename and prefix + m_outobs = (*this)["outobs"].filename(); + m_prefix = (*this)["prefix"].string(); + + // Issue warning if output filename has no .xml suffix + log_string(TERSE, warn_xml_suffix(m_outobs)); + + // Loop over all observation in the container + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[i]); + + // Skip observations that are no CTA observations + if (obs == NULL) { + continue; + } + + // Skip observations that have empty names + if (m_infiles[i].length() == 0) { + continue; + } + + // Set event output file name + std::string outfile = set_outfile_name(m_infiles[i]); + + // Append event extension name + outfile += "["+m_evtname[i]+"]"; + + // Log filename + log_value(NORMAL, "Event list file", outfile); + + // Store output file name in observation + obs->eventfile(outfile); + + // Save event list + save_event_list(obs, m_infiles[i], m_evtname[i], m_gtiname[i], + outfile); + + } // endfor: looped over observations + + // Write observation definition XML file name into logger + log_value(NORMAL, "Obs. definition file", m_outobs); + + // Save observations in XML file + m_obs.save(m_outobs); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save event list into FITS file + * + * @param[in] obs Pointer to CTA observation. + * @param[in] infile Input file name. + * @param[in] evtname Event extension name. + * @param[in] gtiname GTI extension name. + * @param[in] outfile Output file name. + * + * Saves an event list and the corresponding Good Time Intervals into a FITS + * file and copy all others extensions from the input file to the output + * file. + * + * If an extension name is specified in the @p outfile argument, the events + * and eventually also the Good Time Intervals will be extracted from the + * argument and used for writing the events. The format is + * + * [] + * + * where needs to be replaced by the name of the FITS file, + * and by the name of the events + * and Good Time Intervals extensions. For example + * + * myfits.fits[EVENTS1;GTI1] + * + * will write the selected events into the "EVENTS1" extension and the + * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS + * file. If the Good Time Intervals extension name is skipped, e.g. + * + * myfits.fits[EVENTS1] + * + * the original extension name for the Good Time Intervals will be kept. + * Analogously, only the Good Time Intervals extension name can be changed + * by specifying + * + * myfits.fits[;GTI1] + * + * In none of the cases will the original events and Good Time Intervals be + * copied over to the output file. + ***************************************************************************/ +void ctphase::save_event_list(const GCTAObservation* obs, + const std::string& infile, + const std::string& evtname, + const std::string& gtiname, + const std::string& outfile) const +{ + // Save only if we have an event list + if (obs->eventtype() == "EventList") { + + // Set output FITS file event extension names + GFilename outname(outfile); + std::string outevt = evtname; + std::string outgti = gtiname; + if (outname.has_extname()) { + std::vector extnames = + gammalib::split(outname.extname(), ";"); + if (extnames.size() > 0) { + std::string extname = gammalib::strip_whitespace(extnames[0]); + if (!extname.empty()) { + outevt = extname; + } + } + if (extnames.size() > 1) { + std::string extname = gammalib::strip_whitespace(extnames[1]); + if (!extname.empty()) { + outgti = extname; + } + } + } + + // Create output FITS file + GFits outfits; + + // Write observation into FITS file + obs->write(outfits, outevt, outgti); + + // Copy all extensions other than evtname and gtiname extensions + // from the input to the output event list. The evtname and + // gtiname extensions are written by the save method, all others + // that may eventually be present have to be copied over + // explicitly. + GFits infits(infile); + for (int extno = 1; extno < infits.size(); ++extno) { + GFitsHDU* hdu = infits.at(extno); + if (hdu->extname() != evtname && + hdu->extname() != gtiname && + hdu->extname() != outevt && + hdu->extname() != outgti) { + outfits.append(*hdu); + } + } + + // Close input file + infits.close(); + + // Save file to disk and close it (we need both operations) + outfits.saveto(outname.url(), clobber()); + outfits.close(); + + } // endif: observation was unbinned + + // Return + return; +} diff --git a/src/ctphase/ctphase.hpp b/src/ctphase/ctphase.hpp new file mode 100644 index 0000000..dac9f8a --- /dev/null +++ b/src/ctphase/ctphase.hpp @@ -0,0 +1,104 @@ +/*************************************************************************** + * ctphase - Data selection tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2017 by Leonardo Di Venere * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctphase.hpp + * @brief Data selection tool definition + * @author Juergen Knoedlseder + */ + +#ifndef CTPHASE_HPP +#define CTPHASE_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include +#include "GammaLib.hpp" +#include "GCTALib.hpp" +#include "ctobservation.hpp" + +/* __Definitions _________________________________________________________ */ +#define CTPHASE_NAME "ctphase" +#define CTPHASE_VERSION "0.0.1" + + +/***********************************************************************//** + * @class ctphase + * + * @brief Append phase columm to observation files + ***************************************************************************/ +class ctphase : public ctobservation { + +public: + // Constructors and destructors + ctphase(void); + explicit ctphase(const GObservations& obs); + ctphase(int argc, char *argv[]); + ctphase(const ctphase& app); + virtual ~ctphase(void); + + // Operators + ctphase& operator=(const ctphase& app); + + // Methods + void clear(void); + void run(void); + void save(void); + void publish(const std::string& name = ""); + + // TODO: Add Methods for updating model parameters + +protected: + // Protected methods + void init_members(void); + void copy_members(const ctphase& app); + void free_members(void); + void get_parameters(void); + void phase_events(GCTAObservation* obs, + const std::string& filename, + const std::string& evtname, + const std::string& gtiname); + std::string check_infile(const std::string& filename, + const std::string& evtname) const; + std::string set_outfile_name(const std::string& filename) const; + std::string get_gtiname(const std::string& filename, + const std::string& evtname) const; + void save_fits(void); + void save_xml(void); + void save_event_list(const GCTAObservation* obs, + const std::string& infile, + const std::string& evtname, + const std::string& gtiname, + const std::string& outfile) const; + + // User parameters + std::string m_outobs; //!< Output event list or XML file + std::string m_prefix; //!< Prefix for multiple event lists + GModelTemporalPhaseCurve m_phase; //!< Phase information object + GChatter m_chatter; //!< Chattiness + + // Protected members + std::vector m_infiles; //!< Input event filenames + std::vector m_evtname; //!< Event extension names + std::vector m_gtiname; //!< GTI extension names +}; + + +#endif /* CTPHASE_HPP */ diff --git a/src/ctphase/ctphase.par b/src/ctphase/ctphase.par new file mode 100644 index 0000000..9ee197b --- /dev/null +++ b/src/ctphase/ctphase.par @@ -0,0 +1,41 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: ctphase.par # +# Version: 0.0.1 # +# # +# Author: Leonardo Di Venere # +# leonardo.divenere@ba.infn.it # +# INFN # +# # +# Purpose: Parameter file for the ctphase application # +# # +############################################################################# +# +# General parameters +#=================== +inobs, f, a, "events.fits",,, "Input event list or observation definition XML file" +outobs, f, a, "phased_events.fits",,, "Output event list or observation definition XML file" +prefix, s, h, "phased_",,, "Prefix for event lists in observation definition XML file" + +# +# Selection parameters +#===================== +mjd, r, a, 52400,,, "Reference MJD for phase information" +f0, r, a, 1,,, "Frequency at reference MJD" +f1, r, a, 0,,, "First freq. derivative at reference MJD" +f2, r, a, 0,,, "Second freq. derivative at reference MJD" +p0, r, a, 0.1,0.0,1.0, "Phase at reference MJD" + +# +# Standard parameters +#==================== +publish, b, h, no,,, "Publish event list on VO Hub?" +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, "h",,, "Mode of automatic parameters" +logfile, f, h, "ctphase.log",,, "Log filename" diff --git a/src/ctphase/main.cpp b/src/ctphase/main.cpp new file mode 100644 index 0000000..3afb1ab --- /dev/null +++ b/src/ctphase/main.cpp @@ -0,0 +1,77 @@ +/*************************************************************************** + * ctphase - Append phase information to CTA events file * + * ----------------------------------------------------------------------- * + * copyright (C) 2017 by Leonardo Di Venere * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctphase/main.cpp + * @brief Append phase information to CTA events file + * @author Leonardo Di Venere + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "support.hpp" +#include "ctphase.hpp" + + +/***********************************************************************//** + * @brief Main entry point + * + * @param[in] argc Number of arguments + * @param[in] argv Arguments + * + * This is the main entry point of the ctphase application. It allocates a + * ctphase object and executes the application. Any exceptions that occur + * will be catched and corresponding error messages written in the + * application logger and into the standard output. + ***************************************************************************/ +int main (int argc, char *argv[]) +{ + // Initialise return code + int rc = 1; + + // Initialise pointer on application + ctphase* application = NULL; + + // Execute application + try { + + // Create instance of application + application = new ctphase(argc, argv); + + // Execute ctool + rc = execute_ctool(application); + + } + + catch (std::exception &e) { + + // Report exception + report_ctool_failure("ctphase", e.what()); + + } + + // Delete application + delete application; + + // Return + return rc; +} -- libgit2 0.21.4