diff --git a/src/ctskymap/Makefile.am b/src/ctskymap/Makefile.am new file mode 100644 index 0000000..a7a9975 --- /dev/null +++ b/src/ctskymap/Makefile.am @@ -0,0 +1,21 @@ +## Process this file with automake to produce Makefile.in + +# Program name +bin_PROGRAMS = ctskymap + +# Sources +sources = ctskymap.cpp + +# Libtool library +noinst_LTLIBRARIES = libctskymap.la +libctskymap_la_SOURCES = $(sources) + +# Program definition +ctskymap_SOURCES = main.cpp +ctskymap_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@ +ctskymap_CPPFLAGS = @CPPFLAGS@ +ctskymap_LDADD = libctskymap.la + +# Parameter files +pfilesdir = $(prefix)/syspfiles +pfiles_DATA = ctskymap.par diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp new file mode 100644 index 0000000..0137f4d --- /dev/null +++ b/src/ctskymap/ctskymap.cpp @@ -0,0 +1,587 @@ +/*************************************************************************** + * ctskymap - CTA sky mapping tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2010-2011 by Jurgen Knodlseder * + * ----------------------------------------------------------------------- * + * * + * 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 ctskymap.cpp + * @brief CTA sky mapping tool implementation + * @author J. Knodlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctskymap.hpp" +#include "GTools.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_BIN_EVENTS "ctskymap::bin_events(GCTAObservation*)" + +/* __ Debug definitions __________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +ctskymap::ctskymap(void) : GApplication(CTSKYMAP_NAME, CTSKYMAP_VERSION) +{ + // Initialise members + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Observations constructor + * + * This method creates an instance of the class by copying an existing + * observations container. + ***************************************************************************/ +ctskymap::ctskymap(GObservations obs) : GApplication(CTSKYMAP_NAME, CTSKYMAP_VERSION) +{ + // Initialise members + init_members(); + + // Set observations + m_obs = obs; + + // Write header into logger + log_header(); + + // Return + return; +} + + + +/***********************************************************************//** + * @brief Command line constructor + * + * @param[in] argc Number of arguments in command line. + * @param[in] argv Array of command line arguments. + ***************************************************************************/ +ctskymap::ctskymap(int argc, char *argv[]) : + GApplication(CTSKYMAP_NAME, CTSKYMAP_VERSION, argc, argv) +{ + // Initialise members + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] app Application. + ***************************************************************************/ +ctskymap::ctskymap(const ctskymap& app) : GApplication(app) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +ctskymap::~ctskymap(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] app Application. + ***************************************************************************/ +ctskymap& ctskymap::operator= (const ctskymap& app) +{ + // Execute only if object is not identical + if (this != &app) { + + // Copy base class members + this->GApplication::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 instance + ***************************************************************************/ +void ctskymap::clear(void) +{ + // Free members + free_members(); + this->GApplication::free_members(); + + // Initialise members + this->GApplication::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Execute application + * + * This method bins the events data into a counts maps and saves the counts + * map into a FITS file. + ***************************************************************************/ +void ctskymap::execute(void) +{ + // Read ahead output filename so that it gets dumped correctly in the + // parameters log + m_outfile = (*this)["outfile"].filename(); + + // Bin the event data + run(); + + // Save the counts map into FITS file + save(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Bin the event data + ***************************************************************************/ +void ctskymap::run(void) +{ + // Switch screen logging on in debug mode + if (logDebug()) + log.cout(true); + + // Get parameters + get_parameters(); + + // Write parameters into logger + if (logTerse()) { + log_parameters(); + log << std::endl; + } + + // Write observation(s) into logger + if (logTerse()) { + log << std::endl; + if (m_obs.size() > 1) + log.header1("Observations"); + else + log.header1("Observation"); + log << m_obs << std::endl; + } + + // Write header + if (logTerse()) { + log << std::endl; + if (m_obs.size() > 1) + log.header1("Bin observations"); + else + log.header1("Bin observation"); + } + + // 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]); + + // Continue only if observation is a CTA observation + if (obs != NULL) { + + // Write header for observation + if (logTerse()) { + if (obs->name().length() > 1) + log.header3("Observation "+obs->name()); + else + log.header3("Observation"); + } + + // Map events into sky map + map_events(obs); + + } // endif: CTA observation found + + } // endfor: looped over observations + + // Write observation(s) into logger + if (logTerse()) { + log << std::endl; + if (m_obs.size() > 1) + log.header1("Binned observations"); + else + log.header1("Binned observation"); + log << m_obs << std::endl; + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save observation + * + * This method saves the counts map(s) into (a) FITS file(s). + ***************************************************************************/ +void ctskymap::save(void) +{ +/* + // Write header + if (logTerse()) { + log << std::endl; + if (m_obs.size() > 1) + log.header1("Save observations"); + else + log.header1("Save observation"); + } + + // Get output filename + m_outfile = (*this)["outfile"].filename(); + + // Loop over all observations in the container + int file_num = 0; + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(&m_obs[i]); + + // Save only if observation is a CTA observation + if (obs != NULL) { + + // Set filename. If more than one file will be created an + // index "_xxx" will be appended. + std::string filename = m_outfile; + if (file_num > 0) + filename += "_"+str(file_num); + + // Save file + obs->save(filename, clobber()); + + // Increment file number + file_num++; + + } // endif: observation was a CTA observation + + } // endfor: looped over files +*/ + // Return + return; +} + + +/***********************************************************************//** + * @brief Get application parameters + * + * Get all task parameters from parameter file or (if required) by querying + * the user. Most parameters are only required if no observation exists so + * far in the observation container. In this case, a single CTA observation + * will be added to the container, using the definition provided in the + * parameter file. + ***************************************************************************/ +void ctskymap::get_parameters(void) +{ + // If there are no observations in container then add a single CTA + // observation using the parameters from the parameter file + if (m_obs.size() == 0) { + + // Get name of CTA events file + m_evfile = (*this)["evfile"].filename(); + + // Load unbinned CTA observation + GCTAObservation obs; + obs.load_unbinned(m_evfile); + + // Append CTA observation to container + m_obs.append(obs); + + // Use the xref and yref parameters for binning (otherwise the + // pointing direction(s) is/are used) + m_xref = (*this)["xref"].real(); + m_yref = (*this)["yref"].real(); + + } // endif: there was no observation in the container + + // Get remaining parameters + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); + m_proj = (*this)["proj"].string(); + m_coordsys = (*this)["coordsys"].string(); + m_binsz = (*this)["binsz"].real(); + m_nxpix = (*this)["nxpix"].integer(); + m_nypix = (*this)["nypix"].integer(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Map events into a sky map + * + * @param[in] obs CTA observation. + * + * @exception GException::no_list + * No event list found in observation. + * @exception GCTAException::no_pointing + * No valid CTA pointing found. + * + * This method maps the events found in a CTA events list into a sky map. + ***************************************************************************/ +void ctskymap::map_events(GCTAObservation* obs) +{ + // Continue only if observation pointer is valid + if (obs != NULL) { + + // Make sure that the observation holds a CTA event list. If this + // is not the case then throw an exception. + if (dynamic_cast(obs->events()) == NULL) { + throw GException::no_list(G_BIN_EVENTS); + } + + // Setup energy range covered by data + GEnergy emin; + GEnergy emax; + GEbounds ebds; + emin.TeV(m_emin); + emax.TeV(m_emax); + + // Get map centre + double xref; + double yref; + if (m_xref != 9999.0 && m_yref != 9999.0) { + xref = m_xref; + yref = m_yref; + } + else { + + // Get pointer on CTA pointing + const GTime srcTime; // Dummy time + const GCTAPointing *pnt = obs->pointing(srcTime); + if (pnt == NULL) + throw GCTAException::no_pointing(G_BIN_EVENTS); + + // Set reference point to pointing + if (toupper(m_coordsys) == "GAL") { + xref = pnt->dir().l_deg(); + yref = pnt->dir().b_deg(); + } + else { + xref = pnt->dir().ra_deg(); + yref = pnt->dir().dec_deg(); + } + + } // endelse: map centre set to pointing + + // Create skymap + GSkymap map = GSkymap(m_proj, m_coordsys, + xref, yref, m_binsz, m_binsz, + m_nxpix, m_nypix, 1); + + // Initialise binning statistics + int num_outside_map = 0; + int num_outside_erange = 0; + int num_in_map = 0; + + // Fill sky map + GCTAEventList* events = (GCTAEventList*)obs->events(); + for (GCTAEventList::iterator event = events->begin(); event != events->end(); ++event) { + + // Skip if energy is out of range + if (event->energy() < emin || event->energy() > emax) { + num_outside_erange++; + continue; + } + + // Determine sky pixel + GCTAInstDir* inst = (GCTAInstDir*)&(event->dir()); + GSkyDir dir = inst->skydir(); + GSkyPixel pixel = map.dir2xy(dir); + + // Skip if pixel is out of range + if (pixel.x() < -0.5 || pixel.x() > (m_nxpix-0.5) || + pixel.y() < -0.5 || pixel.y() > (m_nypix-0.5)) { + num_outside_map++; + continue; + } + + // Fill event in skymap + map(pixel) += 1.0; + num_in_map++; + + } // endfor: looped over all events + + // Log binning results + if (logTerse()) { + log << std::endl; + log.header1("Mapping"); + log << parformat("Events in list"); + log << obs->events()->size() << std::endl; + log << parformat("Events in map"); + log << num_in_map << std::endl; + log << parformat("Events outside map area"); + log << num_outside_map << std::endl; + log << parformat("Events outside energy range"); + log << num_outside_erange << std::endl; + } + + // Log map + if (logTerse()) { + log << std::endl; + log.header1("Sky map"); + log << map << std::endl; + } + + } // endif: observation was valid + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctskymap::init_members(void) +{ + // Initialise members + m_evfile.clear(); + m_outfile.clear(); + m_proj.clear(); + m_coordsys.clear(); + m_obs.clear(); + m_emin = 0.0; + m_emax = 0.0; + m_xref = 9999.0; // Flags unset (use pointing direction) + m_yref = 9999.0; // Flags unset (use pointing direction) + m_binsz = 0.0; + m_nxpix = 0; + m_nypix = 0; + + // Set logger properties + log.date(true); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctskymap::copy_members(const ctskymap& app) +{ + // Copy attributes + m_evfile = app.m_evfile; + m_outfile = app.m_outfile; + m_proj = app.m_proj; + m_coordsys = app.m_coordsys; + m_obs = app.m_obs; + m_emin = app.m_emin; + m_emax = app.m_emax; + m_xref = app.m_xref; + m_yref = app.m_yref; + m_binsz = app.m_binsz; + m_nxpix = app.m_nxpix; + m_nypix = app.m_nypix; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctskymap::free_members(void) +{ + // Write separator into logger + if (logTerse()) + log << std::endl; + + // Return + return; +} diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp new file mode 100644 index 0000000..f226968 --- /dev/null +++ b/src/ctskymap/ctskymap.hpp @@ -0,0 +1,90 @@ +/*************************************************************************** + * ctskymap - CTA sky mapping tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2010-2011 by Jurgen Knodlseder * + * ----------------------------------------------------------------------- * + * * + * 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 ctskymap.hpp + * @brief CTA sky mapping tool definition + * @author J. Knodlseder + */ + +#ifndef CTSKYMAP_HPP +#define CTSKYMAP_HPP + +/* __ Includes ___________________________________________________________ */ +#include "GammaLib.hpp" +#include "GCTALib.hpp" + +/* __Definitions _________________________________________________________ */ +#define CTSKYMAP_NAME "ctskymap" +#define CTSKYMAP_VERSION "00-01-00" + + +/***********************************************************************//** + * @class ctskymap + * + * @brief CTA sky mapping tool interface defintion + * + * This class creates a sky map from a CTA event list. + ***************************************************************************/ +class ctskymap : public GApplication { +public: + // Constructors and destructors + ctskymap(void); + explicit ctskymap(GObservations obs); + ctskymap(int argc, char *argv[]); + ctskymap(const ctskymap& app); + virtual ~ctskymap(void); + + // Operators + ctskymap& operator= (const ctskymap& app); + + // Methods + void clear(void); + void execute(void); + void run(void); + void save(void); + GObservations& obs(void) { return m_obs; } + void get_parameters(void); + void map_events(GCTAObservation* obs); + +protected: + // Protected methods + void init_members(void); + void copy_members(const ctskymap& app); + void free_members(void); + + // User parameters + std::string m_evfile; //!< Input event list + std::string m_outfile; //!< Output counts map + double m_emin; //!< Minimum energy (TeV) + double m_emax; //!< Maximum energy (TeV) + std::string m_proj; //!< WCS projection + std::string m_coordsys; //!< Coordinate system + double m_xref; //!< Longitude reference coordinate + double m_yref; //!< Latitude reference coordinate + double m_binsz; //!< Pixel size + int m_nxpix; //!< Number of pixels in longitude + int m_nypix; //!< Number of pixels in latitude + + // Protected members + GObservations m_obs; //!< Observation container +}; + +#endif /* CTSKYMAP_HPP */ diff --git a/src/ctskymap/ctskymap.par b/src/ctskymap/ctskymap.par new file mode 100644 index 0000000..9a45d2f --- /dev/null +++ b/src/ctskymap/ctskymap.par @@ -0,0 +1,45 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: ctskymap.par # +# Version: 00-01-00 # +# # +# Author: Jurgen Knodlseder # +# knodlseder@cesr.fr # +# IRAP # +# # +# Purpose: Parameter file for the ctskymap application # +# # +############################################################################# +# +# General parameters +#=================== +evfile, f, a, "test/data/crab_events.fits.gz",,,"Event data file name" +outfile, f, a, "skymap.fits",,, "Output file name" +# +# Event selection parameters +#=========================== +emin, r, a, 0.1,,,"Minimum energy in TeV" +emax, r, a, 100.0,,,"Maximum energy in TeV" +# +# Spatial binning parameters +#=========================== +nxpix, i, a, 200,,, "Size of the X axis in pixels" +nypix, i, a, 200,,, "Size of the Y axis in pixels" +binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" +coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" +xref, r, a, 83.63,,, "First coordinate of image center in degrees (RA or galactic l)" +yref, r, a, 22.01,,, "Second coordinate of image center in degrees (DEC or galactic b)" +axisrot, r, h, 0.0,,, "Rotation angle of image axis, in degrees" +proj, s, a, "CAR",,, "Projection method e.g. AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN" + +# +# Standard parameters +#==================== +chatter, i, h, 1,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, "ql",,, "Mode of automatic parameters" diff --git a/src/ctskymap/main.cpp b/src/ctskymap/main.cpp new file mode 100644 index 0000000..ca7351d --- /dev/null +++ b/src/ctskymap/main.cpp @@ -0,0 +1,79 @@ +/*************************************************************************** + * ctskymap - CTA sky mapping main code * + * ----------------------------------------------------------------------- * + * copyright (C) 2010-2011 by Jurgen Knodlseder * + * ----------------------------------------------------------------------- * + * * + * 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 main.cpp + * @brief CTA sky mapping tool main code + * @author J. Knodlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctskymap.hpp" + + +/***********************************************************************//** + * @brief Main entry point of application + * + * @param[in] argc Number of command line arguments. + * @param[in] argv Command line arguments. + * + * This is the main entry point of the ctskymap application. It allocates a + * ctskymap object and runs the application. + ***************************************************************************/ +int main (int argc, char *argv[]) +{ + // Initialise return code + int rc = 1; + + // Create instance of application + ctskymap application(argc, argv); + + // Run application + try { + // Execute application + application.execute(); + + // Signal success + rc = 0; + } + catch (std::exception &e) { + + // Extract error message + std::string message = e.what(); + std::string signal = "*** ERROR encounterted in the execution of" + " ctskymap. Run aborted ..."; + + // Write error in logger + application.log << signal << std::endl; + application.log << message << std::endl; + + // Write error on standard output + std::cout << signal << std::endl; + std::cout << message << std::endl; + + } // endcatch: catched any application error + + // Return + return rc; +}