From db226c73e7d6915f2f3af074044c4cc978f49f23 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 18 Jul 2014 13:51:36 +0200 Subject: [PATCH] ctbin now always merges event lists into a single counts cube. --- ChangeLog | 3 ++- NEWS | 19 +++++++++++++++++-- pyext/ctbin.i | 7 +++++-- src/ctbin/ctbin.cpp | 334 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------ src/ctbin/ctbin.hpp | 57 ++++++++++++++++++++++++++++++++++++++++++++++++++------- src/ctbin/ctbin.par | 2 +- 6 files changed, 397 insertions(+), 25 deletions(-) diff --git a/ChangeLog b/ChangeLog index b124ed9..c2de956 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2014-07-10 Juergen Knoedlseder +2014-07-18 Juergen Knoedlseder * Version 0.8.0 released ======================== + ctbin now always merges event lists into a single counts cube Add ctexpcube tool Support energy dispersion in ctobssim, ctlike and ctmodel Correct ctobssim for thread safety diff --git a/NEWS b/NEWS index d4b648c..4ccdd38 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00 Author(s) : Juergen Knoedlseder -10 July 2014 +18 July 2014 1. Introduction @@ -40,7 +40,22 @@ when iterating over ctselect (#1149). 6. ctbin - CTA event binning ---------------------------- -None +WARNING: The philosophy behind the ctbin tool changed fundamentally! + +ctbin now always produces a single counts cube on output. If a single +event list is specified on input, the functionnality is basically +unchanged. If multiple event lists are specified, however, events are +now merged into a single counts cube file. No observation definition +XML file is written anymore! + +So far, ctbin still writes pointing information in the counts cube, +although this is not necessarily meaningful if merging event lists +with different pointings. In the long run, pointing information will +be summarized in a pointings definition file, and will be removed +from the counts cube header. The logic is that the pointing definition +is set to the map centre. This is only identical to the pointing in +case of binning a single event list and using the usepnt=yes option +(or specifying the pointing manually as the xref and yref parameters). 7. ctlike - CTA maximum likelihood fitting diff --git a/pyext/ctbin.i b/pyext/ctbin.i index cb2fde8..e215496 100644 --- a/pyext/ctbin.i +++ b/pyext/ctbin.i @@ -1,7 +1,7 @@ /*************************************************************************** * ctbin - CTA data binning tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2013 by Juergen Knoedlseder * + * copyright (C) 2010-2014 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -50,8 +50,11 @@ public: void run(void); void save(void); const GObservations& obs(void) const; + GCTAEventCube cube(void) const; void get_parameters(void); - void bin_events(GCTAObservation* obs); + void init_cube(void); + void fill_cube(GCTAObservation* obs); + void obs_cube(void); }; diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index 6d756f6..df25149 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -34,10 +34,14 @@ /* __ Method name definitions ____________________________________________ */ #define G_BIN_EVENTS "ctbin::bin_events(GCTAObservation*)" +#define G_GET_EBOUNDS "ctbin::get_ebounds()" +#define G_INIT_CUBE "ctbin::init_cube()" +#define G_FILL_CUBE "ctbin::fill_cube(GCTAObservation*)" /* __ Debug definitions __________________________________________________ */ /* __ Coding definitions _________________________________________________ */ +#define G_MERGE_EVENTS /*========================================================================== @@ -270,6 +274,11 @@ void ctbin::run(void) // Initialise observation counter int n_observations = 0; + // Compile option: if cube merging is requested, initialise the cube + #if defined(G_MERGE_EVENTS) + init_cube(); + #endif + // Loop over all observations in the container for (int i = 0; i < m_obs.size(); ++i) { @@ -298,13 +307,25 @@ void ctbin::run(void) // Save event file name (for possible saving) m_infiles[i] = obs->eventfile(); - // Bin events into counts map + // Compile option: if cube merging is requested, fill the cube; + // otherwise just bin the events + #if defined(G_MERGE_EVENTS) + fill_cube(obs); + #else bin_events(obs); + #endif } // endif: CTA observation found } // endfor: looped over observations + // Compile option: if cube merging is requested, set a single cube in + // the observation container + #if defined(G_MERGE_EVENTS) + obs_cube(); + n_observations = 1; + #endif + // 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) { @@ -358,6 +379,12 @@ void ctbin::save(void) } } + // Compile option: if merging was requested, always save into a single + // FITS file + #if defined(G_MERGE_EVENTS) + save_fits(); + #else + /* // Case A: Save counts map(s) and XML metadata information if (m_use_xml) { @@ -371,6 +398,7 @@ void ctbin::save(void) // Always save FITS and an XML summary save_xml(); + #endif // Return return; @@ -423,11 +451,6 @@ void ctbin::get_parameters(void) } - // 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 parameters @@ -436,7 +459,12 @@ void ctbin::get_parameters(void) m_xref = (*this)["xref"].real(); m_yref = (*this)["yref"].real(); } - m_ebinalg = (*this)["ebinalg"].string(); + 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(); + m_ebinalg = (*this)["ebinalg"].string(); // If we have the binning given by a file then read filename if (m_ebinalg == "FILE") { @@ -449,11 +477,6 @@ void ctbin::get_parameters(void) m_emax = (*this)["emax"].real(); m_enumbins = (*this)["enumbins"].integer(); } - 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(); // Optionally read ahead parameters so that they get correctly // dumped into the log file @@ -664,6 +687,206 @@ void ctbin::bin_events(GCTAObservation* obs) } +/***********************************************************************//** + * @brief Initialise counts cube information + * + * @exception GException::invalid_value + * No valid CTA observation found to derive the counts cube + * map centre. + * + * Initialises the skymap, energy boundaries and GTI for a counts cube. + ***************************************************************************/ +void ctbin::init_cube(void) +{ + // Initialse cube information + m_ontime = 0.0; + m_livetime = 0.0; + m_cube.clear(); + m_ebounds.clear(); + m_gti.clear(); + + // Set event cube centre, either from the user parameters or from the + // pointing + double xref = m_xref; + double yref = m_yref; + if (m_usepnt) { + + // Dummy: get pointing from first observation. Ultimately, we want + // to get the pointing from a kind of average + bool found = false; + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* obs = dynamic_cast(m_obs[i]); + if (obs != NULL) { + const GCTAPointing& pnt = obs->pointing(); + if (gammalib::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(); + } + found = true; + break; + } + } + + // Signal if no pointing is found + if (!found) { + std::string msg = "No valid CTA observation has been found in " + "observation list, hence no pointing information " + "could be extracted. Use the \"usepnt=no\" " + "option and specify pointing explicitly when " + "running ctbin."; + throw GException::invalid_value(G_INIT_CUBE, msg); + } + + } // endif: used pointing + + // Create skymap + m_cube = GSkymap(m_proj, m_coordsys, + xref, yref, -m_binsz, m_binsz, + m_nxpix, m_nypix, m_enumbins); + + // Set energy boundaries + get_ebounds(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Fill events into counts cube + * + * @param[in] obs CTA observation. + * + * @exception GException::no_list + * No event list found in observation. + * + * Fills the events from an event list in the counts cube setup by init_cube. + ***************************************************************************/ +void ctbin::fill_cube(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); + } + + // Initialise binning statistics + int num_outside_map = 0; + int num_outside_ebds = 0; + int num_in_map = 0; + + // Fill sky map + GCTAEventList* events = + static_cast(const_cast(obs->events())); + for (int i = 0; i < events->size(); ++i) { + + // Get event + GCTAEventAtom* event = (*events)[i]; + + // Determine sky pixel + GCTAInstDir* inst = (GCTAInstDir*)&(event->dir()); + GSkyDir dir = inst->dir(); + GSkyPixel pixel = m_cube.dir2pix(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; + } + + // Determine energy bin. Skip if we are outside the energy range + int index = m_ebounds.index(event->energy()); + if (index == -1) { + num_outside_ebds++; + continue; + } + + // Fill event in skymap + m_cube(pixel, index) += 1.0; + num_in_map++; + + } // endfor: looped over all events + + // Append GTIs + m_gti.extend(events->gti()); + + // Update ontime and livetime + m_ontime += obs->ontime(); + m_livetime += obs->livetime(); + + // Log filling results + if (logTerse()) { + log << gammalib::parformat("Events in list"); + log << obs->events()->size() << std::endl; + log << gammalib::parformat("Events in cube"); + log << num_in_map << std::endl; + log << gammalib::parformat("Events outside cube area"); + log << num_outside_map << std::endl; + log << gammalib::parformat("Events outside energy bins"); + log << num_outside_ebds << std::endl; + } + + // Log cube + if (logExplicit()) { + log.header1("Counts cube"); + log << m_cube << std::endl; + } + + } // endif: observation was valid + + // Return + return; +} + + +/***********************************************************************//** + * @brief Put cube as single element in observation container. + ***************************************************************************/ +void ctbin::obs_cube(void) +{ + // Clear observation container + m_obs.clear(); + + // Allocate CTA observation. We need this to make sure that all + // attributes are set correctly + GCTAObservation obs; + + // Attach event cube to CTA observation + obs.events(this->cube()); + + // Set map centre as pointing + GSkyPixel pixel(0.5*double(m_cube.nx()), 0.5*double(m_cube.ny())); + GSkyDir centre = m_cube.pix2dir(pixel); + GCTAPointing pointing(centre); + + // Compute deadtime correction + double deadc = (m_ontime > 0.0) ? m_livetime / m_ontime : 0.0; + + // Set CTA observation attributes + obs.pointing(pointing); + obs.obs_id(0); + obs.ra_obj(centre.ra_deg()); //!< Dummy + obs.dec_obj(centre.dec_deg()); //!< Dummy + obs.ontime(m_ontime); + obs.livetime(m_livetime); + obs.deadc(deadc); + + // Append observation to container + m_obs.append(obs); + + // Return + return; +} + + /*========================================================================== = = = Private methods = @@ -698,6 +921,11 @@ void ctbin::init_members(void) m_infiles.clear(); m_use_xml = false; m_read_ahead = false; + m_cube.clear(); + m_ebounds.clear(); + m_gti.clear(); + m_ontime = 0.0; + m_livetime = 0.0; // Set logger properties log.date(true); @@ -737,6 +965,11 @@ void ctbin::copy_members(const ctbin& app) m_infiles = app.m_infiles; m_use_xml = app.m_use_xml; m_read_ahead = app.m_read_ahead; + m_cube = app.m_cube; + m_ebounds = app.m_ebounds; + m_gti = app.m_gti; + m_ontime = app.m_ontime; + m_livetime = app.m_livetime; // Return return; @@ -759,6 +992,83 @@ void ctbin::free_members(void) /***********************************************************************//** + * @brief Get the energy boundaries + * + * @exception GException::invalid_value + * No valid energy boundary extension found. + * + * Get the energy boundaries according to the user parameters. The method + * supports loading of energy boundary information from the EBOUNDS or + * ENERGYBINS extension, or setting energy boundaries using a linear + * or logarithmical spacing. + * + * @todo Ultimately this method should go in a support library as it is + * used by several ctools. + ***************************************************************************/ +void ctbin::get_ebounds(void) +{ + // Initialse energy boundary information + m_ebounds.clear(); + + // Check whether energy binning information should be read from a FITS + // file ... + if (m_ebinalg == "FILE") { + + // Open energy boundary file using the EBOUNDS or ENERGYBINS + // extension. Throw an exception if opening fails. + GFits file(m_ebinfile); + if (file.contains("EBOUNDS")) { + file.close(); + m_ebounds.load(m_ebinfile,"EBOUNDS"); + } + else if (file.contains("ENERGYBINS")) { + file.close(); + m_ebounds.load(m_ebinfile,"ENERGYBINS"); + } + else { + file.close(); + std::string msg = "No extension with name \"EBOUNDS\" or" + " \"ENERGYBINS\" found in FITS file" + " \""+m_ebinfile+"\".\n" + "An \"EBOUNDS\" or \"ENERGYBINS\" extension" + " is required if the parameter \"ebinalg\"" + " is set to \"FILE\"."; + throw GException::invalid_value(G_BIN_EVENTS, msg); + } + + // Set enumbins parameter to number of ebounds + m_enumbins = m_ebounds.size(); + + } // endif: ebinalg was "FILE" + + // ... otherwise use a linear or a logarithmically-spaced energy binning + else { + + // Initialise log mode for ebinning + bool log = true; + + // check if algorithm is linear + if (m_ebinalg == "LIN") { + log = false; + } + + // todo: should we also check if m_ebinalg is "LOG" + // and throw an exception if neither LIN/LOG/FILE + // is given? + + // Setup energy range covered by data + GEnergy emin(m_emin, "TeV"); + GEnergy emax(m_emax, "TeV"); + m_ebounds = GEbounds(m_enumbins, emin, emax, log); + + } //endif: ebinalg was not "FILE" + + // Return + return; +} + + +/***********************************************************************//** * @brief Set output file name. * * @param[in] index Index of input file. diff --git a/src/ctbin/ctbin.hpp b/src/ctbin/ctbin.hpp index d2ce5bf..9e2729f 100644 --- a/src/ctbin/ctbin.hpp +++ b/src/ctbin/ctbin.hpp @@ -33,7 +33,7 @@ /* __Definitions _________________________________________________________ */ #define CTBIN_NAME "ctbin" -#define CTBIN_VERSION "00-04-00" +#define CTBIN_VERSION "00-05-00" /***********************************************************************//** @@ -41,11 +41,24 @@ * * @brief CTA data binning tool interface defintion * - * This class bins CTA event list(s) into a counts map(s). The class can - * operate on predefined observation containers, on individual event list - * FITS files, and on observation definition XML files. Results are stored - * in an observation container that can be written to disk in form of FITS - * files (counts maps) and an updated observation definition file. + * This class bins CTA event list(s) into a single counts cube. The class + * can operate on predefined observation containers, on individual event list + * FITS files, and on observation definition XML files. + * + * If multiple event lists are specified in the observation container or the + * XML definition file, the class will merge these events into a single + * counts cube. + * + * Results are stored in an observation container that can be written to disk + * in form of a single FITS file. + * + * WARNING: Note that the pointing direction of the counts cube will be set + * to the skymap centre used for the counts cube definition. If usepnt=yes + * is used, the pointing direction will be extracted from the first + * observation encountered in the list. Ultimately, pointing direction + * information should not be removed from the counts cube and exposure and + * PSF cubes should be used for response computation. This is however not + * yet fully implemented. ***************************************************************************/ class ctbin : public GApplication { public: @@ -65,20 +78,27 @@ public: void run(void); void save(void); const GObservations& obs(void) const; + GCTAEventCube cube(void) const; void get_parameters(void); - void bin_events(GCTAObservation* obs); + void init_cube(void); + void fill_cube(GCTAObservation* obs); + void obs_cube(void); protected: // Protected methods void init_members(void); void copy_members(const ctbin& app); void free_members(void); + void get_ebounds(void); std::string set_outfile_name(const int index) const; void save_fits(void); void save_xml(void); void save_counts_map(const GCTAObservation* obs, const std::string& outfile) const; + // Obsolete methods (to be reomved later) + void bin_events(GCTAObservation* obs); + // User parameters std::string m_evfile; //!< Input event list or XML file std::string m_outfile; //!< Output counts map or XML file @@ -102,6 +122,11 @@ protected: std::vector m_infiles; //!< Input event filenames bool m_use_xml; //!< Use XML file instead of FITS file bool m_read_ahead; //!< Read ahead parameters + GSkymap m_cube; //!< Event cube + GEbounds m_ebounds; //!< Energy boundaries + GGti m_gti; //!< Good time intervals + double m_ontime; //!< Total ontime + double m_livetime; //!< Total livetime }; /***********************************************************************//** @@ -115,4 +140,22 @@ const GObservations& ctbin::obs(void) const return m_obs; } + +/***********************************************************************//** + * @brief Return event cube + * + * @return Event cube + * + * Returns a CTA event cube built from the ctbin member information. + ***************************************************************************/ +inline +GCTAEventCube ctbin::cube(void) const +{ + // Build event cube + GCTAEventCube cube(m_cube, m_ebounds, m_gti); + + // Return cube + return cube; +} + #endif /* CTBIN_HPP */ diff --git a/src/ctbin/ctbin.par b/src/ctbin/ctbin.par index 6e9dd3c..95b438d 100644 --- a/src/ctbin/ctbin.par +++ b/src/ctbin/ctbin.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctbin.par # -# Version: 00-04-00 # +# Version: 00-05-00 # # # # Author: Jurgen Knodlseder # # Jurgen.Knodlseder@irap.omp.eu # -- libgit2 0.21.4