diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index a6f36f8..9d34c89 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -70,7 +70,7 @@ ctbin::ctbin(void) : ctool(CTBIN_NAME, CTBIN_VERSION) * * param[in] obs Observation container. * - * Constructs a ctbin tool from an observation container. + * Constructs ctbin tool from an observation container. ***************************************************************************/ ctbin::ctbin(const GObservations& obs) : ctool(CTBIN_NAME, CTBIN_VERSION) { @@ -92,7 +92,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. + * Constructs ctbin tool using command line arguments for user parameter + * setting. ***************************************************************************/ ctbin::ctbin(int argc, char *argv[]) : ctool(CTBIN_NAME, CTBIN_VERSION, argc, argv) @@ -149,7 +150,7 @@ ctbin::~ctbin(void) /***********************************************************************//** * @brief Assignment operator * - * @param[in] app Application. + * @param[in] app ctbin tool. * @return ctbin tool. * * Assigns ctbin tool. @@ -240,14 +241,7 @@ void ctbin::run(void) // Write header for observation if (logTerse()) { - std::string header = m_obs[i]->instrument() + " observation"; - if (!m_obs[i]->name().empty()) { - header += " \"" + m_obs[i]->name() + "\""; - } - if (!m_obs[i]->id().empty()) { - header += " (id=" + m_obs[i]->id() +")"; - } - log.header3(header); + log.header3(get_obs_header(m_obs[i])); } // Get CTA observation @@ -328,8 +322,8 @@ void ctbin::save(void) // Log filename if (logTerse()) { - log << "Save counts cube into file \""+m_outcube+"\"."; - log << std::endl; + log << gammalib::parformat("Counts cube file"); + log << m_outcube.url() << std::endl; } // Save cube diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index 68b664e..4486bea 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -60,6 +60,8 @@ const double g_roi_margin = 0.5; //!< Simulation radius margin (degrees) /***********************************************************************//** * @brief Void constructor + * + * Constructs an empty ctobssim tool. ***************************************************************************/ ctobssim::ctobssim(void) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION) { @@ -76,7 +78,7 @@ ctobssim::ctobssim(void) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION) * * @param[in] obs Observation container. * - * Constructs application by initialising from an observation container. + * Constructs ctobssim tool from an observation container. ***************************************************************************/ ctobssim::ctobssim(const GObservations& obs) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION) @@ -98,7 +100,8 @@ ctobssim::ctobssim(const GObservations& obs) : * @param[in] argc Number of arguments in command line. * @param[in] argv Array of command line arguments. * - * Constructs application using command line arguments for parameter setting. + * Constructs ctobssim tool using command line arguments for user parameter + * setting. ***************************************************************************/ ctobssim::ctobssim(int argc, char *argv[]) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv) @@ -115,6 +118,8 @@ ctobssim::ctobssim(int argc, char *argv[]) : * @brief Copy constructor * * @param[in] app Application. + * + * Constructs ctobssim tool from another ctobssim instance. ***************************************************************************/ ctobssim::ctobssim(const ctobssim& app) : ctool(app) { @@ -131,6 +136,8 @@ ctobssim::ctobssim(const ctobssim& app) : ctool(app) /***********************************************************************//** * @brief Destructor + * + * Destructs ctobssim tool. ***************************************************************************/ ctobssim::~ctobssim(void) { @@ -151,8 +158,10 @@ ctobssim::~ctobssim(void) /***********************************************************************//** * @brief Assignment operator * - * @param[in] app Application. - * @return Application. + * @param[in] app ctobssim tool. + * @return ctobssim tool. + * + * Assigns ctobssim tool. ***************************************************************************/ ctobssim& ctobssim::operator=(const ctobssim& app) { @@ -185,7 +194,9 @@ ctobssim& ctobssim::operator=(const ctobssim& app) ==========================================================================*/ /***********************************************************************//** - * @brief Clear application + * @brief Clear ctobssim tool + * + * Clears ctobssim tool. ***************************************************************************/ void ctobssim::clear(void) { @@ -205,10 +216,10 @@ void ctobssim::clear(void) /***********************************************************************//** - * @brief Simulate event data + * @brief Run the ctobssim tool. * - * This method runs the simulation. Results are not saved by this method. - * Invoke "save" to save the results. + * Gets the user parameters, loops over all CTA observations in the + * observation container, and simulate events for each observation. ***************************************************************************/ void ctobssim::run(void) { @@ -286,24 +297,14 @@ void ctobssim::run(void) // 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("Simulate observations"); - } - else { - log.header1("Simulate observation"); - } + log.header1(gammalib::number("Simulate observation", m_obs.size())); } // From here on the code can be parallelized if OpenMP support @@ -334,14 +335,7 @@ void ctobssim::run(void) // Write header for observation if (logTerse()) { - std::string header = m_obs[i]->instrument() + " observation"; - if (m_obs[i]->name().length() > 1) { - header += " \"" + m_obs[i]->name() + "\""; - } - if (m_obs[i]->id().length() > 1) { - header += " (id=" + m_obs[i]->id() +")"; - } - wrklog.header3(header); + wrklog.header3(get_obs_header(m_obs[i])); } // Get pointer on CTA observation @@ -1455,8 +1449,8 @@ void ctobssim::save_fits(void) // Log filename if (logTerse()) { - log << "Save event list into file \""+m_outevents+"\"."; - log << std::endl; + log << gammalib::parformat("Event list file"); + log << m_outevents << std::endl; } // Save observation into FITS file @@ -1483,7 +1477,7 @@ void ctobssim::save_xml(void) { // Get output filename and prefix m_outevents = (*this)["outevents"].filename(); - m_prefix = (*this)["prefix"].string(); + m_prefix = (*this)["prefix"].string(); // Issue warning if output filename has no .xml suffix std::string suffix = gammalib::tolower(m_outevents.substr(m_outevents.length()-4,4)); @@ -1513,15 +1507,16 @@ void ctobssim::save_xml(void) if (obs->events()->size() != 0) { // Set event output file name - std::string outfile = m_prefix + gammalib::str(i) + ".fits"; + std::string outfile = m_prefix + gammalib::str(i) + + ".fits"; // Store output file name in observation obs->eventfile(outfile); // Log filename if (logTerse()) { - log << "Save event list into file \""+outfile+"\"."; - log << std::endl; + log << gammalib::parformat("Event list file"); + log << outfile << std::endl; } // Save observation into FITS file diff --git a/src/ctobssim/ctobssim.par b/src/ctobssim/ctobssim.par index e2a8101..f6bfb2e 100644 --- a/src/ctobssim/ctobssim.par +++ b/src/ctobssim/ctobssim.par @@ -7,7 +7,7 @@ # File: ctobssim.par # # Version: 1.1.0 # # # -# Author: Jurgen Knodlseder # +# Author: Juergen Knoedlseder # # Jurgen.Knodlseder@irap.omp.eu # # IRAP # # # @@ -17,13 +17,13 @@ # # General parameters #=================== -inobs, f, h, "NONE",,, "Input event list or observation definition XML file" -inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file" -caldb, s, a, "prod2",,, "Calibration database" -irf, s, a, "South_0.5h",,, "Instrument response function" +inobs, f, h, NONE,,, "Input event list or observation definition XML file" +inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model XML file" +caldb, s, a, prod2,,, "Calibration database" +irf, s, a, South_0.5h,,, "Instrument response function" edisp, b, h, no,,, "Apply energy dispersion?" -outevents, f, a, "events.fits",,, "Output event data file or observation definition XML file" -prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition XML file" +outevents, f, a, events.fits,,, "Output event data file or observation definition XML file" +prefix, s, h, sim_events_,,, "Prefix for event lists in observation definition XML file" # # Simulation parameters @@ -38,7 +38,7 @@ emin, r, a, 0.1,,, "Lower energy limit (TeV)" emax, r, a, 100.0,,, "Upper energy limit (TeV)" deadc, r, h, 0.95,0,1, "Average deadtime correction factor" maxrate, r, h, 1.0e6,,, "Maximum photon rate" -eslices, i, h, 10,1,100,"Number of energy slices" +eslices, i, h, 10,1,100, "Number of energy slices" # # Standard parameters @@ -46,5 +46,5 @@ eslices, i, h, 10,1,100,"Number of energy slices" 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, "ctobssim.log",,, "Log filename" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, ctobssim.log,,, "Log filename" diff --git a/test/test_ctbin.py b/test/test_ctbin.py index ec8b3a2..8884f0a 100755 --- a/test/test_ctbin.py +++ b/test/test_ctbin.py @@ -22,9 +22,9 @@ import gammalib import ctools -# ========================= # -# Test class for ctbin tool # -# ========================= # +# =========================== # +# Test class for ctbin module # +# =========================== # class Test(gammalib.GPythonTestSuite): """ Test class for ctbin module. @@ -39,7 +39,7 @@ class Test(gammalib.GPythonTestSuite): gammalib.GPythonTestSuite.__init__(self) # Set members - self.events_name = "data/crab_events.fits" + self._events_name = "data/crab_events.fits" # Return return @@ -65,7 +65,7 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctbin bin = ctools.ctbin() - bin["inobs"] = self.events_name + bin["inobs"] = self._events_name bin["outcube"] = "cntmap.fits" bin["ebinalg"] = "LOG" bin["emin"] = 0.1 @@ -111,20 +111,14 @@ class Test(gammalib.GPythonTestSuite): #Please specify a valid parameter name. # Save counts cube - self.test_try("Save counts cube") - try: - bin.save() - self.test_try_success() - except: - msg = "Exception occured in saving counts cube." - self.test_try_failure(msg) + bin.save() # Load counts cube and check content. evt = gammalib.GCTAEventCube("cntmap.fits") self._test_cube(evt, 5542) # Prepare observation container for stacked analysis - cta = gammalib.GCTAObservation(self.events_name) + cta = gammalib.GCTAObservation(self._events_name) obs = gammalib.GObservations() cta.id("0001") obs.append(cta) diff --git a/test/test_ctobssim.py b/test/test_ctobssim.py index 1b36ae2..5816c06 100755 --- a/test/test_ctobssim.py +++ b/test/test_ctobssim.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This scripts performs unit tests for the ctobssim tool. +# This scripts performs unit tests for the ctobssim module. # # Copyright (C) 2014-2016 Juergen Knoedlseder # @@ -20,15 +20,17 @@ # ========================================================================== import gammalib import ctools +from cscripts import obsutils -# ============================ # -# Test class for ctobssim tool # -# ============================ # +# ============================== # +# Test class for ctobssim module # +# ============================== # class Test(gammalib.GPythonTestSuite): """ - Test class for ctobssim tool. + Test class for ctobssim module. """ + # Constructor def __init__(self): """ @@ -38,9 +40,9 @@ class Test(gammalib.GPythonTestSuite): gammalib.GPythonTestSuite.__init__(self) # Set members - self.model_name = "data/crab.xml" - self.caldb = "irf" - self.irf = "cta_dummy_irf" + self._model_name = "data/crab.xml" + self._caldb = "irf" + self._irf = "cta_dummy_irf" # Return return @@ -54,23 +56,23 @@ class Test(gammalib.GPythonTestSuite): self.name("ctobssim") # Append tests - self.append(self.test_functional, "Test ctobssim functionality") - self.append(self.test_container, "Test ctobssim on observation container") + self.append(self._test_ctobssim, "Test ctobssim") + #self.append(self._test_container, "Test ctobssim on observation container") # Return return - # Test ctobssim functionnality - def test_functional(self): + # Test ctobssim + def _test_ctobssim(self): """ - Test ctobssim functionnality. + Test ctobssim. """ # Set-up ctobssim sim = ctools.ctobssim() - sim["inmodel"] = self.model_name + sim["inmodel"] = self._model_name sim["outevents"] = "events.fits" - sim["caldb"] = self.caldb - sim["irf"] = self.irf + sim["caldb"] = self._caldb + sim["irf"] = self._irf sim["ra"] = 83.63 sim["dec"] = 22.01 sim["rad"] = 10.0 @@ -80,154 +82,123 @@ class Test(gammalib.GPythonTestSuite): sim["emax"] = 100.0 # Run tool - self.test_try("Run ctobssim") - try: - sim.run() - self.test_try_success() - #except Exception as e: - # msg = "Exception occured in ctobssim: %s." % (e,) - except: - msg = "Exception occured in ctobssim." - self.test_try_failure(msg) + sim.run() - # Retrieve observation and check content - obs = gammalib.GCTAObservation(sim.obs()[0]) - pnt = obs.pointing() - evt = obs.events() - self.test_assert(obs.instrument() == "CTA", "Observation not a CTA observation") - self.test_value(sim.obs().size(), 1, "There is not a single observation") - self.test_value(obs.ontime(), 1800.0, 1.0e-6, "Ontime is not 1800 sec") - self.test_value(obs.livetime(), 1710.0, 1.0e-6, "Livetime is not 1710 sec") - self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6, "ROI Right Ascension is not 83.63 deg") - self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6, "ROI Declination is not 22.01 deg") - self.test_value(evt.size(), 4119, "Number of events is not 4119") + # Check content of observation + self._test_observation(sim) + self._test_list(sim.obs()[0].events(), 4119) # Save events - self.test_try("Save events") - try: - sim.save() - self.test_try_success() - #except Exception as e: - # msg = "Exception occured in saving events: %s." % (e,) - except: - msg = "Exception occured in saving events." - self.test_try_failure(msg) + sim.save() + + # Load counts cube and check content. + evt = gammalib.GCTAEventList("events.fits") + self._test_list(evt, 4119) - # Test ctobssim on observation container - def test_container(self): - """ - Test ctobssim on observation container. - """ # Set-up observation container - obs = self.set_obs(4) + pnts = [{'ra': 83.63, 'dec': 21.01}, + {'ra': 84.63, 'dec': 22.01}, + {'ra': 83.63, 'dec': 23.01}, + {'ra': 82.63, 'dec': 22.01}] + obs = obsutils.set_obs_list(pnts, caldb=self._caldb, irf=self._irf) - # Set-up ctobssim + # Set-up ctobssim from observation container sim = ctools.ctobssim(obs) sim["outevents"] = "sim_events.xml" - sim["inmodel"] = self.model_name + sim["inmodel"] = self._model_name + # Run tool - self.test_try("Run ctobssim") - try: - sim.run() - self.test_try_success() - #except Exception as e: - # msg = "Exception occured in ctobssim: %s." % (e,) - except: - msg = "Exception occured in ctobssim." - self.test_try_failure(msg) + sim.run() # Retrieve observation and check content - self.test_value(sim.obs().size(), 4, "There are not 4 observations") + self._test_observation(sim, nobs=4, pnts=pnts) + self._test_list(sim.obs()[0].events(), 4063) + self._test_list(sim.obs()[1].events(), 4059) + self._test_list(sim.obs()[2].events(), 4010) + self._test_list(sim.obs()[3].events(), 4138) # Save events - self.test_try("Save events") + sim.save() + + # Load events + obs = gammalib.GObservations("sim_events.xml") + + # Retrieve observation and check content + self._test_list(obs[0].events(), 4063) + self._test_list(obs[1].events(), 4059) + self._test_list(obs[2].events(), 4010) + self._test_list(obs[3].events(), 4138) + + # Set-up ctobssim with invalid event file + sim = ctools.ctobssim() + sim["inmodel"] = "model_file_that_does_not_exist.xml" + sim["outevents"] = "events.fits" + sim["caldb"] = self._caldb + sim["irf"] = self._irf + sim["ra"] = 83.63 + sim["dec"] = 22.01 + sim["rad"] = 10.0 + sim["tmin"] = 0.0 + sim["tmax"] = 1800.0 + sim["emin"] = 0.1 + sim["emax"] = 100.0 + + # Run ctbin tool + self.test_try("Run ctobssim with invalid model file") try: - sim.save() - self.test_try_success() - #except Exception as e: - # msg = "Exception occured in saving events: %s." % (e,) + sim.run() + self.test_try_failure() except: - msg = "Exception occured in saving events." - self.test_try_failure(msg) + self.test_try_success() + + # Return + return - # Setup observation container - def set_obs(self, number): + # Check observation + def _test_observation(self, ctobssim, nobs=1, + pnts=[{'ra': 83.63, 'dec': 22.01}]): """ - Setup observation container. + Test content of an observation. + + Args: + ctobssim: ctobssim instance + + Kwargs: + nobs: Number of observations + pnts: List of pointing dictionaries """ - # Initialise empty observation container - obs = gammalib.GObservations() - - # Initialise first time and identifier - tstart = 0.0 - offset = -float(number)/2.0 - - # Loop over number of observations - for i in range(number): - id = "%6.6d" % i - obs_cta = self.set_one_obs(id, offset) - offset += 1.0 - tstart += 1830.0 - obs.append(obs_cta) - - # Return observation container - return obs - - # Setup one observation - def set_one_obs(self, id, offset, \ - tstart=0.0, duration=1800.0, deadc=0.95, \ - emin=0.1, emax=100.0, rad=5.0, \ - irf="cta_dummy_irf", caldb="dummy"): + # Test observation container + self.test_value(ctobssim.obs().size(), nobs, + "There is one observation") + for i in range(ctobssim.obs().size()): + obs = gammalib.GCTAObservation(ctobssim.obs()[i]) + pnt = obs.pointing() + 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(), pnts[i]['ra'], 1.0e-6, + "Pointing Right Ascension is "+str(pnts[i]['ra'])+" deg") + self.test_value(pnt.dir().dec_deg(), pnts[i]['dec'], 1.0e-6, + "Pointing Declination is "+str(pnts[i]['dec'])+" deg") + + # Return + return + + # Check event list + def _test_list(self, list, nevents): """ - Setup one observation for test purposes. + Test content of event list." + + Args: + list: Event list + nevents: Expected number of events. """ - # Allocate CTA observation - obs_cta = gammalib.GCTAObservation() - - # Set calibration database - db = gammalib.GCaldb() - if (gammalib.dir_exists(self.caldb)): - db.rootdir(self.caldb) - else: - db.open("cta", self.caldb) - - # Set pointing direction - pntdir = gammalib.GSkyDir() - pntdir.radec_deg(83.63, 22.01+offset) - pnt = gammalib.GCTAPointing() - pnt.dir(pntdir) - obs_cta.pointing(pnt) - - # Set ROI - roi = gammalib.GCTARoi() - instdir = gammalib.GCTAInstDir() - instdir.dir(pntdir) - roi.centre(instdir) - roi.radius(rad) - - # Set GTI - gti = gammalib.GGti() - gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration)) - - # Set energy boundaries - ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \ - gammalib.GEnergy(emax, "TeV")) - - # Allocate event list - events = gammalib.GCTAEventList() - events.roi(roi) - events.gti(gti) - events.ebounds(ebounds) - obs_cta.events(events) - - # Set instrument response - obs_cta.response(self.irf, db) - - # Set ontime, livetime, and deadtime correction factor - obs_cta.ontime(duration) - obs_cta.livetime(duration*deadc) - obs_cta.deadc(deadc) - obs_cta.id(id) - - # Return CTA observation - return obs_cta + # Test event list + self.test_value(list.size(), nevents, str(nevents)+" elements") + self.test_value(list.number(), nevents, str(nevents)+" events") + + # Return + return