From 0847394c57541fad7444abdfc8dc62c6c11cd418 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 15 Mar 2016 00:15:01 +0100 Subject: [PATCH] Enhance ctbin Python module test. --- cscripts/obsutils.py | 92 ++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------ src/ctbin/ctbin.cpp | 136 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------- src/ctbin/ctbin.par | 22 +++++++++++----------- test/test_ctbin.py | 107 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------- 4 files changed, 228 insertions(+), 129 deletions(-) diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py index fc07d24..e176335 100755 --- a/cscripts/obsutils.py +++ b/cscripts/obsutils.py @@ -1,8 +1,7 @@ # ========================================================================== -# This script provides a number of functions that are useful for handling -# CTA observations. +# CTA observation handling support functions. # -# Copyright (C) 2011-2015 Juergen Knoedlseder +# Copyright (C) 2011-2016 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 @@ -31,21 +30,28 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, """ Simulate events for all observations in the container. - Parameters: - obs - Observation container - Keywords: - log - Create log file(s) - debug - Create console dump? - edisp - Apply energy dispersion? - seed - Seed value for simulations (default: 0) - nbins - Number of energy bins (default: 0=unbinned) - binsz - Pixel size for binned simulation (deg/pixel) - npix - Number of pixels in X and Y for binned simulation + Args: + obs: Observation container without events + + Kwargs: + log: Create log file(s) + debug: Create console dump? + chatter: Chatter level + edisp: Apply energy dispersion? + seed: Seed value for simulations (default: 0) + nbins: Number of energy bins (default: 0=unbinned) + binsz: Pixel size for binned simulation (deg/pixel) + npix: Number of pixels in X and Y for binned simulation + proj: Projection for binned simulation + coord: Coordinate system for binned simulation + + Returns: + Observation container, filled with simulated events. """ # Allocate ctobssim application and set parameters sim = ctools.ctobssim(obs) - sim["seed"] = seed + sim["seed"] = seed sim["edisp"] = edisp # Optionally open the log file @@ -85,16 +91,16 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, # Allocate ctbin application and set parameters bin = ctools.ctbin(sim.obs()) - bin["ebinalg"] = "LOG" - bin["emin"] = emin - bin["emax"] = emax + bin["ebinalg"] = "LOG" + bin["emin"] = emin + bin["emax"] = emax bin["enumbins"] = nbins - bin["usepnt"] = True # Use pointing for map centre - bin["nxpix"] = npix - bin["nypix"] = npix - bin["binsz"] = binsz + bin["usepnt"] = True # Use pointing for map centre + bin["nxpix"] = npix + bin["nypix"] = npix + bin["binsz"] = binsz bin["coordsys"] = coord - bin["proj"] = proj + bin["proj"] = proj # Optionally open the log file if log: @@ -350,28 +356,30 @@ def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ # ======================= # def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ emin=0.1, emax=100.0, rad=5.0, \ - irf="South_50h", caldb="prod2", id="000000", instrument="CTA"): + irf="South_50h", caldb="prod2", id="000000"): """ - Returns a single CTA observation containing an empty CTA event list. - By looping over this function you can add CTA observations to the - observation container. - - Parameters: - pntdir - Pointing direction [GSkyDir] - Keywords: - tstart - Start time [seconds] (default: 0.0) - duration - Duration of observation [seconds] (default: 1800.0) - deadc - Deadtime correction factor (default: 0.95) - emin - Minimum event energy [TeV] (default: 0.1) - emax - Maximum event energy [TeV] (default: 100.0) - rad - ROI radius used for analysis [deg] (default: 5.0) - irf - Instrument response function (default: cta_dummy_irf) - caldb - Calibration database path (default: "dummy") - id - Run identifier (default: "000000") - instrument - Intrument (default: "CTA") + Set a single CTA observation. + + The function sets a single CTA observation containing an empty CTA + event list. By looping over this function you can add CTA observations + to the observation container. + + Args: + pntdir - Pointing direction [GSkyDir] + + Kwargs: + tstart: Start time (seconds) (default: 0.0) + duration: Duration of observation (seconds) (default: 1800.0) + deadc: Deadtime correction factor (default: 0.95) + emin: Minimum event energy (TeV) (default: 0.1) + emax: Maximum event energy (TeV) (default: 100.0) + rad: ROI radius used for analysis (deg) (default: 5.0) + irf: Instrument response function (default: "South_50h") + caldb: Calibration database path (default: "prod2") + id: Run identifier (default: "000000") """ # Allocate CTA observation - obs_cta = gammalib.GCTAObservation(instrument) + obs_cta = gammalib.GCTAObservation() # Set calibration database db = gammalib.GCaldb() @@ -397,7 +405,7 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration)) # Set energy boundaries - ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \ + ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), gammalib.GEnergy(emax, "TeV")) # Allocate event list diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index b169956..a6f36f8 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -52,6 +52,8 @@ const GEnergy g_energy_margin(1.0e-12, "TeV"); /***********************************************************************//** * @brief Void constructor + * + * Constructs an empty ctbin tool. ***************************************************************************/ ctbin::ctbin(void) : ctool(CTBIN_NAME, CTBIN_VERSION) { @@ -68,8 +70,7 @@ ctbin::ctbin(void) : ctool(CTBIN_NAME, CTBIN_VERSION) * * param[in] obs Observation container. * - * This method creates an instance of the class by copying an existing - * observations container. + * Constructs a ctbin tool from an observation container. ***************************************************************************/ ctbin::ctbin(const GObservations& obs) : ctool(CTBIN_NAME, CTBIN_VERSION) { @@ -90,6 +91,8 @@ ctbin::ctbin(const GObservations& obs) : ctool(CTBIN_NAME, CTBIN_VERSION) * * @param[in] argc Number of arguments in command line. * @param[in] argv Array of command line arguments. + * + * Constructs ctbin tool from command line arguments. ***************************************************************************/ ctbin::ctbin(int argc, char *argv[]) : ctool(CTBIN_NAME, CTBIN_VERSION, argc, argv) @@ -106,6 +109,8 @@ ctbin::ctbin(int argc, char *argv[]) : * @brief Copy constructor * * @param[in] app Application. + * + * Constructs ctbin tool from another ctbin instance. ***************************************************************************/ ctbin::ctbin(const ctbin& app) : ctool(app) { @@ -122,6 +127,8 @@ ctbin::ctbin(const ctbin& app) : ctool(app) /***********************************************************************//** * @brief Destructor + * + * Destructs ctbin tool. ***************************************************************************/ ctbin::~ctbin(void) { @@ -143,6 +150,9 @@ ctbin::~ctbin(void) * @brief Assignment operator * * @param[in] app Application. + * @return ctbin tool. + * + * Assigns ctbin tool. ***************************************************************************/ ctbin& ctbin::operator=(const ctbin& app) { @@ -175,7 +185,9 @@ ctbin& ctbin::operator=(const ctbin& app) ==========================================================================*/ /***********************************************************************//** - * @brief Clear instance + * @brief Clear ctbin tool + * + * Clears ctbin tool. ***************************************************************************/ void ctbin::clear(void) { @@ -195,12 +207,10 @@ void ctbin::clear(void) /***********************************************************************//** - * @brief Bin the event data + * @brief Run the ctbin tool * - * This method loops over all observations found in the observation conatiner - * and bins all events from the event list(s) into counts map(s). Note that - * each event list is binned in a separate counts map, hence no summing of - * events is done. + * Gets the user parameters and loops over all CTA observations in the + * observation container to bin the events into a single counts cube. ***************************************************************************/ void ctbin::run(void) { @@ -212,33 +222,17 @@ void ctbin::run(void) // Get task 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.header1(gammalib::number("Observation", m_obs.size())); 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"); - } + log.header1(gammalib::number("Bin observation", m_obs.size())); } // Loop over all observations in the container @@ -247,10 +241,10 @@ void ctbin::run(void) // Write header for observation if (logTerse()) { std::string header = m_obs[i]->instrument() + " observation"; - if (m_obs[i]->name().length() > 1) { + if (!m_obs[i]->name().empty()) { header += " \"" + m_obs[i]->name() + "\""; } - if (m_obs[i]->id().length() > 1) { + if (!m_obs[i]->id().empty()) { header += " (id=" + m_obs[i]->id() +")"; } log.header3(header); @@ -293,12 +287,7 @@ void ctbin::run(void) // 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.header1(gammalib::number("Binned observation", m_obs.size())); log << m_obs << std::endl; } @@ -496,6 +485,12 @@ void ctbin::get_parameters(void) m_outcube = (*this)["outcube"].filename(); } + // Write parameters into logger + if (logTerse()) { + log_parameters(); + log << std::endl; + } + // Return return; } @@ -521,8 +516,8 @@ void ctbin::fill_cube(GCTAObservation* obs) const GCTAEventList* events = dynamic_cast(obs->events()); if (events == NULL) { std::string msg = "CTA Observation does not contain an event " - "list. Event list information is needed to " - "fill the counts map."; + "list. An event list is needed to fill the " + "counts cube."; throw GException::invalid_value(G_FILL_CUBE, msg); } @@ -534,7 +529,7 @@ void ctbin::fill_cube(GCTAObservation* obs) // Check for RoI sanity if (!roi.is_valid()) { - std::string msg = "No RoI information found in input observation " + std::string msg = "No valid RoI found in input observation " "\""+obs->name()+"\". Run ctselect to specify " "an RoI for this observation before running " "ctbin."; @@ -624,10 +619,12 @@ void ctbin::fill_cube(GCTAObservation* obs) * @brief Create output observation container. * * Creates an output observation container that combines all input CTA - * observation into a single cube-style observation. All non CTA observations - * present in the observation container are kept. The method furthermore - * conserves any response information in case that a single CTA observation - * is provided. This supports the original binned analysis. + * observation into a single stacked observation. All non-CTA observations + * and all binned CTA observations that were present in the observation + * container are append to the observation container so that they can + * be used by other tools. The method furthermore conserves any response + * information in case that a single CTA observation is provided to support + * binned analysis. ***************************************************************************/ void ctbin::obs_cube(void) { @@ -641,26 +638,37 @@ void ctbin::obs_cube(void) GCTAObservation* obs = dynamic_cast(m_obs[0]); if (obs != NULL) { - // Only change the event type if we had an unbinned observation + // Change the event type if we had an unbinned observation if (obs->eventtype() == "EventList") { + // Assign cube to the observation obs->events(this->cube()); - obs->eventfile(""); - } // endif: observation was unbinned - else { // input observation was binned and skipped - // Create new and empty cube - GCTAEventCube cube = GCTAEventCube(m_cube, m_ebounds, obs->gti()); + } + + // ... otherwise the input observation was binned and hence + // skipped. In that case we simply append an empty counts + // cube + else { + + // Create empty counts cube + GCTAEventCube cube(m_cube, m_ebounds, obs->gti()); - // Assign empty cube to have a new binned observation + // Assign empty cube to observation obs->events(cube); + } + // Reset file name + obs->eventfile(""); + } // endif: obervation was valid } // endif: we only had one observation in the container - // ... otherwise put a single CTA observation in container + // ... otherwise put a single CTA observation in container and + // append all observations that have not been used in the binning + // to the container else { // Allocate observation container @@ -672,10 +680,25 @@ void ctbin::obs_cube(void) // 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 average pointing direction for all CTA event lists + double ra = 0.0; + double dec = 0.0; + double number = 0.0; + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* cta = dynamic_cast(m_obs[i]); + if ((cta != NULL) && (cta->eventtype() == "EventList")) { + ra += cta->pointing().dir().ra(); + dec += cta->pointing().dir().dec(); + number += 1.0; + } + } + if (number > 0.0) { + ra /= number; + dec /= number; + } + GSkyDir dir; + dir.radec(ra, dec); + GCTAPointing pointing(dir); // Compute deadtime correction double deadc = (m_ontime > 0.0) ? m_livetime / m_ontime : 0.0; @@ -683,8 +706,8 @@ void ctbin::obs_cube(void) // 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.ra_obj(dir.ra_deg()); //!< Dummy + obs.dec_obj(dir.dec_deg()); //!< Dummy obs.ontime(m_ontime); obs.livetime(m_livetime); obs.deadc(deadc); @@ -701,6 +724,9 @@ void ctbin::obs_cube(void) if (obs == NULL) { container.append(*m_obs[i]); } + else if (obs->eventtype() != "EventList") { + container.append(*m_obs[i]); + } } // Set observation container diff --git a/src/ctbin/ctbin.par b/src/ctbin/ctbin.par index 8ce64dc..a55b087 100644 --- a/src/ctbin/ctbin.par +++ b/src/ctbin/ctbin.par @@ -7,27 +7,27 @@ # File: ctbin.par # # Version: 1.1.0 # # # -# Author: Jurgen Knodlseder # +# Author: Juergen Knoedlseder # # Jurgen.Knodlseder@irap.omp.eu # # IRAP # # # -# Purpose: Parameter file for the ctbin application # +# Purpose: Parameter file for the ctbin tool # # # ############################################################################# # # General parameters #=================== -inobs, f, a, "events.fits",,, "Input event list or observation definition XML file" -outcube, f, a, "cntcube.fits",,, "Output counts cube file" +inobs, f, a, events.fits,,, "Input event list or observation definition XML file" +outcube, f, a, cntcube.fits,,, "Output counts cube file" # # Energy binning parameters #========================== -ebinalg, s, a, "LOG",FILE|LIN|LOG,, "Algorithm for defining energy bins" +ebinalg, s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining energy bins" emin, r, a, 0.1,,, "Start value for first energy bin in TeV" emax, r, a, 100.0,,, "Stop value for last energy bin in TeV" -enumbins, i, a, 20,,, "Number of energy bins" -ebinfile, f, a, "NONE",,, "Name of the file containing the energy bin definition" +enumbins, i, a, 20,1,200, "Number of energy bins" +ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition" # # Spatial binning parameters @@ -36,10 +36,10 @@ usepnt, b, h, no,,, "Use pointing instead of xref/yref 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)" +coordsys, s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" -proj, s, a, "CAR",AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" +proj, s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" # # Standard parameters @@ -48,5 +48,5 @@ publish, b, h, no,,, "Publish counts cube 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, "ql",,, "Mode of automatic parameters" -logfile, f, h, "ctbin.log",,, "Log filename" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, ctbin.log,,, "Log filename" diff --git a/test/test_ctbin.py b/test/test_ctbin.py index f17d62d..ec8b3a2 100755 --- a/test/test_ctbin.py +++ b/test/test_ctbin.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This scripts performs unit tests for the ctbin tool. +# This scripts performs unit tests for the ctbin module. # # Copyright (C) 2014-2016 Juergen Knoedlseder # @@ -27,8 +27,9 @@ import ctools # ========================= # class Test(gammalib.GPythonTestSuite): """ - Test class for ctbin tool. + Test class for ctbin module. """ + # Constructor def __init__(self): """ @@ -78,14 +79,8 @@ class Test(gammalib.GPythonTestSuite): bin["xref"] = 83.63 bin["yref"] = 22.01 - # Run tool - self.test_try("Run ctbin") - try: - bin.run() - self.test_try_success() - except: - msg = "Exception occured in ctbin." - self.test_try_failure(msg) + # Run ctbin tool + bin.run() # Check content of observation and cube self._test_observation(bin, 5542) @@ -98,7 +93,7 @@ class Test(gammalib.GPythonTestSuite): self._test_observation(cpy_bin, 5542) self._test_cube(cpy_bin.cube(), 5542) - # Run copy of ctbin again + # Run copy of ctbin tool again cpy_bin.run() # Check content of observation and cube. We expect now an empty @@ -128,46 +123,116 @@ class Test(gammalib.GPythonTestSuite): evt = gammalib.GCTAEventCube("cntmap.fits") self._test_cube(evt, 5542) + # Prepare observation container for stacked analysis + cta = gammalib.GCTAObservation(self.events_name) + obs = gammalib.GObservations() + cta.id("0001") + obs.append(cta) + cta.id("0002") + obs.append(cta) + cta.id("0003") + obs.append(cta) + + # Set-up ctbin using an observation container + bin = ctools.ctbin(obs) + bin["outcube"] = "cntmap.fits" + bin["ebinalg"] = "LOG" + bin["emin"] = 0.1 + bin["emax"] = 100.0 + bin["enumbins"] = 20 + bin["nxpix"] = 200 + bin["nypix"] = 200 + bin["binsz"] = 0.02 + bin["coordsys"] = "CEL" + bin["proj"] = "CAR" + bin["xref"] = 83.63 + bin["yref"] = 22.01 + + # Run ctbin tool + bin.run() + + # Check content of observation and cube (need multiplier=3 since + # three identical observations have been appended) + self._test_observation(bin, 5542, multiplier=3) + self._test_cube(bin.cube(), 5542, multiplier=3) + + # Set-up ctbin with invalid event file + bin = ctools.ctbin() + bin["inobs"] = "event_file_that_does_not_exist.fits" + bin["outcube"] = "cntmap.fits" + bin["ebinalg"] = "LOG" + bin["emin"] = 0.1 + bin["emax"] = 100.0 + bin["enumbins"] = 20 + bin["nxpix"] = 200 + bin["nypix"] = 200 + bin["binsz"] = 0.02 + bin["coordsys"] = "CEL" + bin["proj"] = "CAR" + bin["xref"] = 83.63 + bin["yref"] = 22.01 + + # Run ctbin tool + self.test_try("Run ctbin with invalid event file") + try: + bin.run() + self.test_try_failure() + except: + self.test_try_success() + # Return return # Check observation - def _test_observation(self, ctbin, nevents): + def _test_observation(self, ctbin, nevents, multiplier=1): """ Test content of an observation. Args: ctbin: ctbin instance nevents: Expected number of events + + Kwargs: + multiplier: Observation multiplier """ # Test observation container obs = gammalib.GCTAObservation(ctbin.obs()[0]) pnt = obs.pointing() - self.test_value(ctbin.obs().size(), 1, "There is one observation") - self.test_assert(obs.instrument() == "CTA", "Observation is CTA observation") - self.test_value(obs.ontime(), 1800.0, 1.0e-6, "Ontime is 1800 sec") - self.test_value(obs.livetime(), 1710.0, 1.0e-6, "Livetime is 1710 sec") - self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6, "Pointing Right Ascension is 83.63 deg") - self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6, "Pointing Declination is 22.01 deg") + self.test_value(ctbin.obs().size(), 1, + "There is one observation") + self.test_assert(obs.instrument() == "CTA", + "Observation is CTA observation") + self.test_value(obs.ontime(), 1800.0*multiplier, 1.0e-6, + "Ontime is 1800 sec") + self.test_value(obs.livetime(), 1710.0*multiplier, 1.0e-6, + "Livetime is 1710 sec") + self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6, + "Pointing Right Ascension is 83.63 deg") + self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6, + "Pointing Declination is 22.01 deg") # Test event cube - self._test_cube(obs.events(), nevents) + self._test_cube(obs.events(), nevents, multiplier=multiplier) # Return return # Check event cube - def _test_cube(self, cube, nevents): + def _test_cube(self, cube, nevents, multiplier=1): """ Test content of event cube." Args: cube: Event cube nevents: Expected number of events. + + Kwargs: + multiplier: Event number multiplier """ # Test event cube self.test_value(cube.size(), 800000, "800000 event bins") - self.test_value(cube.number(), nevents, str(nevents)+" events") + self.test_value(cube.number(), nevents*multiplier, + str(nevents*multiplier)+" events") # Return return -- libgit2 0.21.4