diff --git a/pyext/ctmapcube.i b/pyext/ctmapcube.i index a94abb5..27b06b6 100644 --- a/pyext/ctmapcube.i +++ b/pyext/ctmapcube.i @@ -44,11 +44,12 @@ public: virtual ~ctmapcube(void); // Methods - void clear(void); - void run(void); - void save(void); - void publish(const std::string& name = ""); - void models(const GModels& models); + void clear(void); + void run(void); + void save(void); + void publish(const std::string& name = ""); + const GModelSpatialDiffuseCube& mapcube(void) const; + void models(const GModels& models); }; diff --git a/pyext/ctool.i b/pyext/ctool.i index ddb7888..5120da1 100644 --- a/pyext/ctool.i +++ b/pyext/ctool.i @@ -84,26 +84,31 @@ public: virtual void execute(void); // Make methods private in Python by prepending an underscore - %rename(_read_ahead) read_ahead; - %rename(_time_reference) time_reference; - %rename(_get_observations) get_observations; - %rename(_setup_observations) setup_observations; - %rename(_create_ebounds) create_ebounds; - %rename(_create_map) create_map; - %rename(_create_cube) create_cube; - %rename(_create_cta_obs) create_cta_obs; - %rename(_require_inobs) require_inobs; - %rename(_require_inobs_nolist) require_inobs_nolist; - %rename(_require_inobs_nocube) require_inobs_nocube; - %rename(_log_observations) log_observations; - %rename(_log_models) log_models; - %rename(_set_response) set_response; - %rename(_set_obs_response) set_obs_response; - %rename(_set_obs_bounds) set_obs_bounds; - %rename(_get_mean_pointing) get_mean_pointing; - %rename(_get_current_rss) get_current_rss; - %rename(_get_obs_header) get_obs_header; - %rename(_warn_too_few_energies) warn_too_few_energies; + %rename(_read_ahead) read_ahead; + %rename(_time_reference) time_reference; + %rename(_get_observations) get_observations; + %rename(_setup_observations) setup_observations; + %rename(_create_ebounds) create_ebounds; + %rename(_create_map) create_map; + %rename(_create_cube) create_cube; + %rename(_create_cta_obs) create_cta_obs; + %rename(_require_inobs) require_inobs; + %rename(_require_inobs_nolist) require_inobs_nolist; + %rename(_require_inobs_nocube) require_inobs_nocube; + %rename(_log_observations) log_observations; + %rename(_log_models) log_models; + %rename(_set_response) set_response; + %rename(_set_edisp) set_edisp; + %rename(_restore_edisp) restore_edisp; + %rename(_set_obs_response) set_obs_response; + %rename(_set_obs_bounds) set_obs_bounds; + %rename(_get_mean_pointing) get_mean_pointing; + %rename(_get_current_rss) get_current_rss; + %rename(_get_obs_header) get_obs_header; + %rename(_insert_energy_boundaries) insert_energy_boundaries; + %rename(_cube_layer_usage) cube_layer_usage; + %rename(_warn_too_few_energies) warn_too_few_energies; + %rename(_warn_xml_suffix) warn_xml_suffix; // Protected methods const bool& read_ahead(void) const; @@ -131,16 +136,23 @@ public: const std::string& what = "Model"); // Protected support methods - void set_response(GObservations& obs); - void set_obs_response(GCTAObservation* obs); - void set_obs_bounds(GObservations& obs); - GSkyDir get_mean_pointing(const GObservations& obs); - size_t get_current_rss(void); - std::string get_obs_header(const GObservation* obs); - GEnergies insert_energy_boundaries(const GEnergies& energies, - const GCTAObservation& obs); + void set_response(GObservations& obs); + std::vector set_edisp(GObservations& obs, const bool& edisp) const; + void restore_edisp(GObservations& obs, + const std::vector& edisp) const; + void set_obs_response(GCTAObservation* obs); + void set_obs_bounds(GObservations& obs); + GSkyDir get_mean_pointing(const GObservations& obs); + size_t get_current_rss(void); + std::string get_obs_header(const GObservation* obs); + GEnergies insert_energy_boundaries(const GEnergies& energies, + const GCTAObservation& obs); + std::vector cube_layer_usage(const GEbounds& cube_ebounds, + const GEbounds& list_ebounds) const; + // Protected warning strings std::string warn_too_few_energies(const GEnergies& energies) const; + std::string warn_xml_suffix(const GFilename& filename) const; }; diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index d2bb22e..ec799d9 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -254,7 +254,7 @@ void ctbin::run(void) // Skip observation if we have a binned observation if (obs->eventtype() == "CountsCube") { - std::string msg = " Skipping binned "+m_obs[i]->instrument()+ + std::string msg = " Skipping binned "+obs->instrument()+ " observation"; log_string(NORMAL, msg); continue; @@ -518,7 +518,7 @@ void ctbin::fill_cube(GCTAObservation* obs) } // Get counts cube usage flags - std::vector usage = cube_layer_usage(events->ebounds()); + std::vector usage = cube_layer_usage(m_ebounds, events->ebounds()); // Initialise binning statistics int num_outside_roi = 0; @@ -624,7 +624,7 @@ void ctbin::set_weights(GCTAObservation* obs) } // Get counts cube usage flags - std::vector usage = cube_layer_usage(events->ebounds()); + std::vector usage = cube_layer_usage(m_ebounds, events->ebounds()); // Loop over all pixels in counts cube for (int pixel = 0; pixel < m_counts.npix(); ++pixel) { @@ -778,42 +778,3 @@ void ctbin::obs_cube(void) // Return return; } - - -/***********************************************************************//** - * @brief Determine the counts cube layer usage. - * - * @param[in] ebounds Energy boundaries of the event list. - * @return Vector of usage flags. - * - * Determines a vector of counts cube layer usage flags that signal whether - * an event should actually be filled in the counts cube or not. This makes - * sure that no partially filled bins will exist in the counts cube. - ***************************************************************************/ -std::vector ctbin::cube_layer_usage(const GEbounds& ebounds) const -{ - // Set energy margin - const GEnergy energy_margin(0.01, "GeV"); - - // Initialise usage vector - std::vector usage(m_ebounds.size(), true); - - // Loop over all energy bins of the cube - for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) { - - // If the counts cube energy bin is fully contained within the - // energy boundaries of the event list the signal usage of this - // counts cube bin. Partially overlapping energy bins are signaled - // for non-usage, which avoids having partially filled bins in the - // counts cube. Some margin is applied that effectively reduces the - // width of the counts cube energy bin. This should cope with any - // rounding errors that come from reading and writing the energy - // boundary information to a FITS file. - usage[iebin] = ebounds.contains(m_ebounds.emin(iebin) + energy_margin, - m_ebounds.emax(iebin) - energy_margin); - - } // endfor: looped over all energy bins - - // Return usage - return usage; -} diff --git a/src/ctbin/ctbin.hpp b/src/ctbin/ctbin.hpp index 138f6c0..73930d1 100644 --- a/src/ctbin/ctbin.hpp +++ b/src/ctbin/ctbin.hpp @@ -85,14 +85,13 @@ public: protected: // Protected methods - void init_members(void); - void copy_members(const ctbin& app); - void free_members(void); - void get_parameters(void); - void fill_cube(GCTAObservation* obs); - void set_weights(GCTAObservation* obs); - void obs_cube(void); - std::vector cube_layer_usage(const GEbounds& ebounds) const; + void init_members(void); + void copy_members(const ctbin& app); + void free_members(void); + void get_parameters(void); + void fill_cube(GCTAObservation* obs); + void set_weights(GCTAObservation* obs); + void obs_cube(void); // User parameters GFilename m_outcube; //!< Output counts map file name diff --git a/src/ctmapcube/ctmapcube.cpp b/src/ctmapcube/ctmapcube.cpp index 3832f36..e1c304f 100644 --- a/src/ctmapcube/ctmapcube.cpp +++ b/src/ctmapcube/ctmapcube.cpp @@ -192,16 +192,16 @@ void ctmapcube::run(void) // Get task parameters get_parameters(); - // Write models into logger - if (logTerse()) { - log << std::endl; - log.header1(gammalib::number("Model", m_models.size())); - log << m_models << std::endl; - } + // Write input model container into logger + log_models(NORMAL, m_models, "Input model"); // Generate map cube create_cube(); + // Write result cube into logger + log_header1(NORMAL, "Map cube"); + log_string(NORMAL, m_cube.print(m_chatter)); + // Optionally publish map cube if (m_publish) { publish(); @@ -220,27 +220,27 @@ void ctmapcube::run(void) void ctmapcube::save(void) { // Write header - if (logTerse()) { - log << std::endl; - log.header1("Save map cube"); - } + log_header1(TERSE, "Save map cube"); // Get map cube filename m_outcube = (*this)["outcube"].filename(); - // Save only if filename is non-empty - if (!m_outcube.is_empty()) { - - // Log filename - if (logTerse()) { - log << gammalib::parformat("Map cube file"); - log << m_outcube.url() << std::endl; - } + // Determine whether the cube is empty + bool cube_is_empty = ((m_cube.cube().nx() == 0) || + (m_cube.cube().ny() == 0) || + (m_cube.cube().nmaps() == 0)); - // Save map cube into FITS file + // Save map cube if filename and the map cube are not empty + if (!m_outcube.is_empty() && !cube_is_empty) { m_cube.save(m_outcube, clobber()); + } + // Write into logger what has been done + std::string fname = (m_outcube.is_empty()) ? "NONE" : m_outcube.url(); + if (cube_is_empty) { + fname.append(" (cube is empty, no file created)"); } + log_value(NORMAL, "Map cube file", fname); // Return return; @@ -256,11 +256,8 @@ void ctmapcube::save(void) ***************************************************************************/ void ctmapcube::publish(const std::string& name) { - // Write header - if (logTerse()) { - log << std::endl; - log.header1("Publish map cube"); - } + // Write header into logger + log_header1(TERSE, "Publish map cube"); // Set default name is user name is empty std::string user_name(name); @@ -268,11 +265,8 @@ void ctmapcube::publish(const std::string& name) user_name = CTMAPCUBE_NAME; } - // Log filename - if (logTerse()) { - log << gammalib::parformat("Map cube"); - log << user_name << std::endl; - } + // Write map cube name into logger + log_value(NORMAL, "Map cube name", user_name); // Publish map cube m_cube.cube().publish(user_name); @@ -297,6 +291,7 @@ void ctmapcube::init_members(void) m_outcube.clear(); m_ptsrcsig = 1.0; m_publish = false; + m_chatter = static_cast(2); // Initialise protected members m_models.clear(); @@ -318,6 +313,7 @@ void ctmapcube::copy_members(const ctmapcube& app) m_outcube = app.m_outcube; m_ptsrcsig = app.m_ptsrcsig; m_publish = app.m_publish; + m_chatter = app.m_chatter; // Copy protected members m_models = app.m_models; @@ -372,6 +368,7 @@ void ctmapcube::get_parameters(void) // Get remaining hidden parameters m_ptsrcsig = (*this)["ptsrcsig"].real(); m_publish = (*this)["publish"].boolean(); + m_chatter = static_cast((*this)["chatter"].integer()); // Read optionally output cube filenames if (read_ahead()) { @@ -462,11 +459,8 @@ GEnergies ctmapcube::create_energies(void) ***************************************************************************/ void ctmapcube::create_cube(void) { - // Write header - if (logTerse()) { - log << std::endl; - log.header1("Generate map cube"); - } + // Write header into logger + log_header1(TERSE, "Generate map cube"); // Loop over all models for (int i = 0; i < m_models.size(); ++i) { @@ -478,10 +472,7 @@ void ctmapcube::create_cube(void) if (model == NULL) { // Signal that model is skipped - if (logNormal()) { - log << gammalib::parformat("Skip model"); - log << m_models[i]->name() << std::endl; - } + log_value(NORMAL, "Skip model", m_models[i]->name()); // Continue continue; @@ -489,10 +480,7 @@ void ctmapcube::create_cube(void) } // Signal that model is used - if (logNormal()) { - log << gammalib::parformat("Use model"); - log << m_models[i]->name() << std::endl; - } + log_value(NORMAL, "Use model", m_models[i]->name()); // Add sky model if (dynamic_cast(model->spatial()) != NULL) { @@ -504,13 +492,6 @@ void ctmapcube::create_cube(void) } // endfor: looped over energies - // Log cube - if (logTerse()) { - log << std::endl; - log.header1("Map cube"); - log << m_cube << std::endl; - } - // Return return; } diff --git a/src/ctmapcube/ctmapcube.hpp b/src/ctmapcube/ctmapcube.hpp index 37d72d2..6bb5b39 100644 --- a/src/ctmapcube/ctmapcube.hpp +++ b/src/ctmapcube/ctmapcube.hpp @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTMAPCUBE_NAME "ctmapcube" -#define CTMAPCUBE_VERSION "1.0.0" +#define CTMAPCUBE_VERSION "1.2.0" /***********************************************************************//** @@ -55,11 +55,12 @@ public: ctmapcube& operator=(const ctmapcube& app); // Methods - void clear(void); - void run(void); - void save(void); - void publish(const std::string& name = ""); - void models(const GModels& models); + void clear(void); + void run(void); + void save(void); + void publish(const std::string& name = ""); + const GModelSpatialDiffuseCube& mapcube(void) const; + void models(const GModels& models); protected: // Protected methods @@ -79,6 +80,7 @@ protected: GFilename m_outcube; //!< Output map cube filename double m_ptsrcsig; //!< Point source sigma (arcmin) bool m_publish; //!< Publish map cube? + GChatter m_chatter; //!< Chattiness // Protected members GModels m_models; //!< Model container @@ -87,6 +89,18 @@ protected: /***********************************************************************//** + * @brief Return map cube + * + * @return Map cube + ***************************************************************************/ +inline +const GModelSpatialDiffuseCube& ctmapcube::mapcube(void) const +{ + return (m_cube); +} + + +/***********************************************************************//** * @brief Set models * * @param[in] models Model container. diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp index 6a79c86..af053e2 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -34,6 +34,7 @@ /* __ Method name definitions ____________________________________________ */ #define G_GET_PARAMETERS "ctmodel::get_parameters()" +#define G_FILL_CUBE "ctmodel::fill_cube(GCTAObservation*)" /* __ Debug definitions __________________________________________________ */ @@ -217,68 +218,33 @@ void ctmodel::run(void) // Get task parameters get_parameters(); - // Set energy dispersion flag for all CTA observations and save old + // Set energy dispersion flags of all CTA observations and save old // values in save_edisp vector - std::vector save_edisp; - save_edisp.assign(m_obs.size(), false); - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - save_edisp[i] = obs->response()->apply_edisp(); - obs->response()->apply_edisp(m_apply_edisp); - } - } + std::vector save_edisp = set_edisp(m_obs, m_apply_edisp); - // 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 input observation container into logger + log_observations(NORMAL, m_obs, "Input observation"); - // Write models into logger - if (logTerse()) { - log << std::endl; - log.header1("Models"); - log << m_obs.models() << std::endl; - } + // Write input model container into logger + log_models(NORMAL, m_obs.models(), "Input model"); // Write header - if (logTerse()) { - log << std::endl; - log.header1("Generate model cube"); - } + log_header1(TERSE, "Generate model cube"); // Loop over all observations in the container for (int i = 0; i < m_obs.size(); ++i) { - // 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() +")"; - } - log.header3(header); - } + // Write header for the current observation + log_header3(TERSE, get_obs_header(m_obs[i])); // Get CTA observation GCTAObservation* obs = dynamic_cast(m_obs[i]); // Skip observation if it's not CTA if (obs == NULL) { - if (logTerse()) { - log << " Skipping "; - log << m_obs[i]->instrument(); - log << " observation" << std::endl; - } + std::string msg = " Skipping "+m_obs[i]->instrument()+ + " observation"; + log_string(NORMAL, msg); continue; } @@ -291,11 +257,9 @@ void ctmodel::run(void) // Skip observation if we have a binned observation if (obs->eventtype() == "CountsCube") { - if (logTerse()) { - log << " Skipping binned "; - log << obs->instrument(); - log << " observation" << std::endl; - } + std::string msg = " Skipping binned "+m_obs[i]->instrument()+ + " observation"; + log_string(NORMAL, msg); continue; } @@ -309,20 +273,12 @@ void ctmodel::run(void) } // endfor: looped over observations - // Log cube - if (logTerse()) { - log << std::endl; - log.header1("Model cube"); - log << m_cube << std::endl; - } + // Write model cube into header + log_header1(NORMAL, "Model cube"); + log_string(NORMAL, m_cube.print()); - // Restore energy dispersion flag for all CTA observations - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - obs->response()->apply_edisp(save_edisp[i]); - } - } + // Restore energy dispersion flags of all CTA observations + restore_edisp(m_obs, save_edisp); // Optionally publish model cube if (m_publish) { @@ -343,27 +299,22 @@ void ctmodel::run(void) void ctmodel::save(void) { // Write header - if (logTerse()) { - log << std::endl; - log.header1("Save cube"); - } + log_header1(TERSE, "Save model cube"); // Get model cube filename m_outcube = (*this)["outcube"].filename(); - // Save only if filename is non-empty - if (!m_outcube.is_empty()) { - - // Log filename - if (logTerse()) { - log << gammalib::parformat("Model cube file"); - log << m_outcube.url() << std::endl; - } - - // Save model cube into FITS file + // Save only if filename is non-empty and cube has some size + if (!m_outcube.is_empty() && m_cube.size() > 0) { m_cube.save(m_outcube, clobber()); + } + // Write into logger what has been done + std::string fname = (m_outcube.is_empty()) ? "NONE" : m_outcube.url(); + if (m_cube.size() == 0) { + fname.append(" (cube is empty, no file created)"); } + log_value(NORMAL, "Model cube file", fname); // Return return; @@ -377,11 +328,8 @@ void ctmodel::save(void) ***************************************************************************/ void ctmodel::publish(const std::string& name) { - // Write header - if (logTerse()) { - log << std::endl; - log.header1("Publish model cube"); - } + // Write header into logger + log_header1(TERSE, "Publish model cube"); // Set default name is user name is empty std::string user_name(name); @@ -390,10 +338,7 @@ void ctmodel::publish(const std::string& name) } // Log filename - if (logTerse()) { - log << gammalib::parformat("Publish model cube"); - log << user_name << std::endl; - } + log_value(NORMAL, "Model cube name", user_name); // Publish model cube m_cube.counts().publish(user_name); @@ -443,6 +388,7 @@ void ctmodel::init_members(void) m_outcube.clear(); m_apply_edisp = false; m_publish = false; + m_chatter = static_cast(2); // Initialise protected members m_obs.clear(); @@ -468,6 +414,7 @@ void ctmodel::copy_members(const ctmodel& app) m_outcube = app.m_outcube; m_apply_edisp = app.m_apply_edisp; m_publish = app.m_publish; + m_chatter = app.m_chatter; // Copy protected members m_obs = app.m_obs; @@ -525,7 +472,8 @@ void ctmodel::get_parameters(void) if (obs->eventtype() == "CountsCube") { // Set cube from binned observation - GCTAEventCube* evtcube = dynamic_cast(const_cast(obs->events())); + GCTAEventCube* evtcube = dynamic_cast + (const_cast(obs->events())); cube(*evtcube); @@ -593,6 +541,7 @@ void ctmodel::get_parameters(void) // Get remaining parameters m_publish = (*this)["publish"].boolean(); + m_chatter = static_cast((*this)["chatter"].integer()); // Read optionally output cube filenames if (read_ahead()) { @@ -776,97 +725,85 @@ void ctmodel::get_obs(void) ***************************************************************************/ void ctmodel::fill_cube(const GCTAObservation* obs) { - // Continue only if observation pointer is valid - if (obs != NULL) { - - // Get GTI and energy boundaries references for observation - const GGti& gti = obs->events()->gti(); - const GEbounds& obs_ebounds = obs->ebounds(); - - // Get cube energy boundaries - const GEbounds& cube_ebounds = m_cube.ebounds(); + // Get references to GTI and energy boundaries for the event list + const GGti& gti = obs->events()->gti(); + const GEbounds& obs_ebounds = obs->ebounds(); - // Initialise empty, invalid RoI - GCTARoi roi; + // Get cube energy boundaries + const GEbounds& cube_ebounds = m_cube.ebounds(); - // Retrieve RoI in case we have an unbinned observation - if (obs->eventtype() == "EventList") { - roi = obs->roi(); - } + // Get counts cube usage flags + std::vector usage = cube_layer_usage(cube_ebounds, obs_ebounds); - // Initialise statistics - double sum = 0.0; - int num_outside_ebds = 0; - int num_outside_roi = 0; + // Initialise empty, invalid RoI + GCTARoi roi; - // Setup cube GTIs for this observation - m_cube.gti(obs->events()->gti()); + // Retrieve RoI in case we have an unbinned observation + if (obs->eventtype() == "EventList") { + roi = obs->roi(); + } - // Loop over all cube bins - for (int i = 0; i < m_cube.size(); ++i) { + // Initialise statistics + double sum = 0.0; + int num_outside_ebds = 0; + int num_outside_roi = 0; - // Get cube bin - GCTAEventBin* bin = m_cube[i]; + // Setup cube GTIs for this observation + m_cube.gti(obs->events()->gti()); - // Skip bin if it is outside the energy range of the observation - int index = cube_ebounds.index(bin->energy()); - if (index == -1 || - !obs_ebounds.contains(cube_ebounds.emin(index)+g_energy_margin, - cube_ebounds.emax(index)-g_energy_margin)) { - num_outside_ebds++; - continue; - } + // Loop over all cube bins + for (int i = 0; i < m_cube.size(); ++i) { - // Check if RoI is valid, i.e. check if we have an unbinned - // observation - if (roi.is_valid()) { + // Get cube bin + GCTAEventBin* bin = m_cube[i]; - // Skip bin if it is outside the RoI of the observation - if (!roi.contains(*bin)) { - num_outside_roi++; - continue; - } + // Determine counts cube energy bin + int iebin = cube_ebounds.index(bin->energy()); - } // endif: RoI was not valid + // Skip bin if the corresponding counts cube energy bin is not fully + // contained in the event list energy range. This avoids having + // partially filled bins. + if (!usage[iebin]) { + num_outside_ebds++; + continue; + } - // Get actual bin value - double value = bin->counts(); - - // Compute model value for cube bin - double model = m_obs.models().eval(*bin, *obs) * bin->size(); + // If RoI is valid then skip bin if it is outside the RoI of the + // observation + if (roi.is_valid() && !roi.contains(*bin)) { + num_outside_roi++; + continue; + } - // Add model to actual value - value += model; - sum += model; + // Get actual bin value + double value = bin->counts(); + + // Compute model value for cube bin + double model = m_obs.models().eval(*bin, *obs) * bin->size(); - // Store value - bin->counts(value); + // Add model to actual value + value += model; + sum += model; - } // endfor: looped over all cube bins + // Store value + bin->counts(value); - // Append GTIs of observation to list of GTIs - m_gti.extend(gti); + } // endfor: looped over all cube bins - // Update GTIs - m_cube.gti(m_gti); + // Append GTIs of observation to list of GTIs + m_gti.extend(gti); - // Log results - if (logTerse()) { - log << gammalib::parformat("Model events in cube"); - log << sum << std::endl; - log << gammalib::parformat("Bins outside energy range"); - log << num_outside_ebds << std::endl; - log << gammalib::parformat("Bins outside RoI"); - log << num_outside_roi << std::endl; - } + // Update GTIs + m_cube.gti(m_gti); - // Log cube - if (logExplicit()) { - log.header2("Model cube"); - log << m_cube << std::endl; - } + // Log filling results + log_value(NORMAL, "Model events in cube", sum); + log_value(NORMAL, "Bins outside energy range", num_outside_ebds); + log_value(NORMAL, "Bins outside RoI", num_outside_roi); - } // endif: observation was valid + // Write model cube into header + log_header2(EXPLICIT, "Model cube"); + log_string(EXPLICIT, m_cube.print(m_chatter)); // Return return; diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp index 5fd44c5..fe2cebc 100644 --- a/src/ctmodel/ctmodel.hpp +++ b/src/ctmodel/ctmodel.hpp @@ -88,6 +88,7 @@ protected: GFilename m_outcube; //!< Output model cube bool m_apply_edisp; //!< Apply energy dispersion? bool m_publish; //!< Publish model cube? + GChatter m_chatter; //!< Chattiness // Protected members GObservations m_obs; //!< Observation container diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index d6302ce..b592ed5 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -242,16 +242,15 @@ void ctobssim::run(void) m_save_and_dispose = true; } - // Determine the number of valid CTA observations, set energy dispersion - // flag for all CTA observations and save old values in save_edisp vector - int n_observations = 0; - std::vector save_edisp; - save_edisp.assign(m_obs.size(), false); + // Set energy dispersion flags of all CTA observations and save old + // values in save_edisp vector + std::vector save_edisp = set_edisp(m_obs, m_apply_edisp); + + // Determine the number of valid CTA observations + int n_observations = 0; for (int i = 0; i < m_obs.size(); ++i) { GCTAObservation* obs = dynamic_cast(m_obs[i]); if (obs != NULL) { - save_edisp[i] = obs->response()->apply_edisp(); - obs->response()->apply_edisp(m_apply_edisp); n_observations++; } } @@ -262,48 +261,29 @@ void ctobssim::run(void) m_use_xml = true; } - // Write execution mode into logger - if (logTerse()) { - log << std::endl; - log.header1("Execution mode"); - log << gammalib::parformat("Event list management"); - if (m_save_and_dispose) { - log << "Save and dispose (reduces memory needs)" << std::endl; - } - else { - log << "Keep events in memory" << std::endl; - } - log << gammalib::parformat("Output format"); - if (m_use_xml) { - log << "Write Observation Definition XML file" << std::endl; - } - else { - log << "Write single event list FITS file" << std::endl; - } - } + // Write execution model into logger + log_header1(NORMAL, "Execution mode"); + std::string mode = (m_save_and_dispose) + ? "Save and dispose (reduces memory needs)" + : "Keep events in memory"; + std::string xml = (m_use_xml) + ? "Write Observation Definition XML file" + : "Write single event list FITS file"; + log_value(NORMAL, "Event list management", mode); + log_value(NORMAL, "Output format", xml); // Write seed values into logger - if (logTerse()) { - log << std::endl; - log.header1("Seed values"); - for (int i = 0; i < m_rans.size(); ++i) { - log << gammalib::parformat("Seed "+gammalib::str(i)); - log << gammalib::str(m_rans[i].seed()) << std::endl; - } + log_header1(NORMAL, "Seed values"); + for (int i = 0; i < m_rans.size(); ++i) { + log_value(NORMAL, "Seed "+gammalib::str(i), + gammalib::str(m_rans[i].seed())); } - // Write observation(s) into logger - if (logTerse()) { - log << std::endl; - log.header1(gammalib::number("Observation", m_obs.size())); - log << m_obs << std::endl; - } + // Write input observation container into logger + log_observations(NORMAL, m_obs, "Input observation"); // Write header - if (logTerse()) { - log << std::endl; - log.header1(gammalib::number("Simulate observation", m_obs.size())); - } + log_header1(TERSE, gammalib::number("Simulate observation", m_obs.size())); // From here on the code can be parallelized if OpenMP support // is enabled. The code in the following block corresponds to the @@ -433,13 +413,8 @@ void ctobssim::run(void) } // end pragma omp parallel - // Restore energy dispersion flag for all CTA observations - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - obs->response()->apply_edisp(save_edisp[i]); - } - } + // Restore energy dispersion flags of all CTA observations + restore_edisp(m_obs, save_edisp); // Return return; @@ -467,16 +442,8 @@ void ctobssim::run(void) ***************************************************************************/ void ctobssim::save(void) { - // Write header - if (logTerse()) { - log << std::endl; - if (m_obs.size() > 1) { - log.header1("Save observations"); - } - else { - log.header1("Save observation"); - } - } + // Write header into logger + log_header1(TERSE, gammalib::number("Save observation", m_obs.size())); // Case A: Save event file(s) and XML metadata information if (m_use_xml) { @@ -1460,10 +1427,7 @@ void ctobssim::save_fits(void) GCTAObservation* obs = dynamic_cast(m_obs[0]); // Log filename - if (logTerse()) { - log << gammalib::parformat("Event list file"); - log << m_outevents << std::endl; - } + log_value(NORMAL, "Event list file", m_outevents); // Save observation into FITS file obs->save(m_outevents, clobber()); @@ -1492,16 +1456,7 @@ void ctobssim::save_xml(void) 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)); - if (suffix != ".xml") { - log << "*** WARNING: Name of observation definition output file \""+ - m_outevents+"\"" << std::endl; - log << "*** WARNING: does not terminate with \".xml\"." << std::endl; - log << "*** WARNING: This is not an error, but might be misleading." - " It is recommended" << std::endl; - log << "*** WARNING: to use the suffix \".xml\" for observation" - " definition files." << std::endl; - } + log_string(TERSE, warn_xml_suffix(m_outevents)); // Save only if event lists have not yet been saved and disposed if (!m_save_and_dispose) { @@ -1526,10 +1481,7 @@ void ctobssim::save_xml(void) obs->eventfile(outfile); // Log filename - if (logTerse()) { - log << gammalib::parformat("Event list file"); - log << outfile << std::endl; - } + log_value(NORMAL, "Event list file", outfile); // Save observation into FITS file obs->save(outfile, clobber()); diff --git a/src/ctselect/ctselect.cpp b/src/ctselect/ctselect.cpp index 29c3cf4..5323eb8 100644 --- a/src/ctselect/ctselect.cpp +++ b/src/ctselect/ctselect.cpp @@ -224,18 +224,11 @@ void ctselect::run(void) // Get parameters get_parameters(); - // Write observation(s) into logger - if (logTerse()) { - log << std::endl; - log.header1("Observations before selection"); - log << m_obs << std::endl; - } + // Write input observation container into logger + log_observations(NORMAL, m_obs, "Input observation"); - // Write header - if (logTerse()) { - log << std::endl; - log.header1("Event selection"); - } + // Write header into logger + log_header1(TERSE, "Event selection"); // Initialise counters int n_observations = 0; @@ -243,17 +236,8 @@ void ctselect::run(void) // Loop over all observation in the container for (int i = 0; i < m_obs.size(); ++i) { - // 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() +")"; - } - log.header3(header); - } + // 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 @@ -266,21 +250,17 @@ void ctselect::run(void) // Skip observation if it's not CTA if (obs == NULL) { - if (logTerse()) { - log << " Skipping "; - log << m_obs[i]->instrument(); - log << " observation" << std::endl; - } + 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") { - if (logTerse()) { - log << " Skipping binned "; - log << obs->instrument(); - log << " observation" << std::endl; - } + std::string msg = " Skipping binned "+m_obs[i]->instrument()+ + " observation"; + log_string(NORMAL, msg); continue; } @@ -297,23 +277,15 @@ void ctselect::run(void) } m_gtiname[i] = get_gtiname(fname.url(), m_evtname[i]); - // Log input file information - if (logTerse()) { - log << gammalib::parformat("Input filename"); - log << m_infiles[i] << std::endl; - log << gammalib::parformat("Event extension name"); - log << m_evtname[i] << std::endl; - log << gammalib::parformat("GTI extension name"); - log << m_gtiname[i] << std::endl; - } + // 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) { - if (logTerse()) { - log << " Warning: No events in event file \""; - log << m_infiles[i] << "\". Event selection skipped."; - log << std::endl; - } + log_string(NORMAL, " Warning: No events in event file \""+ + m_infiles[i]+"\". Event selection skipped."); continue; } @@ -377,11 +349,7 @@ void ctselect::run(void) } // Write observation(s) into logger - if (logTerse()) { - log << std::endl; - log.header1("Observations after selection"); - log << m_obs << std::endl; - } + log_observations(NORMAL, m_obs, "Output observation"); // Return return; @@ -409,16 +377,8 @@ void ctselect::run(void) ***************************************************************************/ void ctselect::save(void) { - // Write header - if (logTerse()) { - log << std::endl; - if (m_obs.size() > 1) { - log.header1("Save event lists"); - } - else { - log.header1("Save event list"); - } - } + // 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) { @@ -459,6 +419,7 @@ void ctselect::init_members(void) m_emax = 0.0; m_expr.clear(); m_usethres.clear(); + m_chatter = static_cast(2); // Initialise protected members m_obs.clear(); @@ -496,6 +457,7 @@ void ctselect::copy_members(const ctselect& app) m_emax = app.m_emax; m_expr = app.m_expr; m_usethres = app.m_usethres; + m_chatter = app.m_chatter; // Copy protected members m_obs = app.m_obs; @@ -610,6 +572,7 @@ void ctselect::get_parameters(void) // Get other User parameters m_expr = (*this)["expr"].string(); m_usethres = (*this)["usethres"].string(); + m_chatter = static_cast((*this)["chatter"].integer()); // Optionally read ahead parameters so that they get correctly // dumped into the log file @@ -726,10 +689,8 @@ void ctselect::select_events(GCTAObservation* obs, sprintf(cmin, "%.8f", tmin); sprintf(cmax, "%.8f", tmax); selection = "TIME >= "+std::string(cmin)+" && TIME <= "+std::string(cmax); - if (logTerse()) { - log << gammalib::parformat("Time range"); - log << tmin << " - " << tmax << " s" << std::endl; - } + log_value(NORMAL, "Time range", + gammalib::str(tmin)+" - "+gammalib::str(tmax)+" s"); } // endif: made time selection @@ -742,7 +703,7 @@ void ctselect::select_events(GCTAObservation* obs, } // Log the requested energy selection - if (logTerse()) { + if (logNormal()) { log << gammalib::parformat("Selected energy range"); if (emin > 0.0 && emax > 0.0) { if (emin >= emax) { @@ -804,23 +765,20 @@ void ctselect::select_events(GCTAObservation* obs, } // Log the requested ROI - if (logTerse()) { - log << gammalib::parformat("Requested ROI"); - log << "Centre(RA,DEC)=(" << ra << ", " << dec << ") deg, "; - log << "Radius=" << m_rad << " deg" << std::endl; - } + log_value(NORMAL, "Requested RoI", + "Centre(RA,DEC)=("+gammalib::str(ra)+", "+ + gammalib::str(dec)+") deg, Radius="+gammalib::str(m_rad)+ + " deg"); // If we have already an ROI then make sure that the selected // ROI overlaps with the existing ROI double roi_radius = list->roi().radius(); if (roi_radius > 0.0) { GSkyDir roi_centre = list->roi().centre().dir(); - if (logTerse()) { - log << gammalib::parformat("ROI of data"); - log << "Centre(RA,DEC)=(" << roi_centre.ra_deg() << ", "; - log << roi_centre.dec_deg() << ") deg, "; - log << "Radius=" << roi_radius << " deg" << std::endl; - } + log_value(NORMAL, "RoI of data", + "Centre(RA,DEC)=("+gammalib::str(roi_centre.ra_deg())+ + ", "+gammalib::str(roi_centre.dec_deg())+") deg, Radius="+ + gammalib::str(roi_radius)+" deg"); GSkyDir centre; centre.radec_deg(ra, dec); double distance = centre.dist_deg(roi_centre); @@ -830,8 +788,8 @@ void ctselect::select_events(GCTAObservation* obs, } // Log the selected ROI - if (logTerse()) { - log << gammalib::parformat("Selected ROI"); + if (logNormal()) { + log << gammalib::parformat("Selected RoI"); if (rad <= 0.0) { log << "None. There is no overlap between existing "; log << "and requested ROI." << std::endl; @@ -869,10 +827,7 @@ void ctselect::select_events(GCTAObservation* obs, } // endif: made expression selection // Dump cfitsio selection string - if (logTerse()) { - log << gammalib::parformat("cfitsio selection"); - log << selection << std::endl; - } + log_value(NORMAL, "cfitsio selection", selection); // Build input filename including selection expression std::string expression = filename; @@ -881,20 +836,14 @@ void ctselect::select_events(GCTAObservation* obs, } // Dump FITS filename including selection expression - if (logTerse()) { - log << gammalib::parformat("FITS filename"); - log << expression << std::endl; - } + log_value(NORMAL, "FITS filename", expression); // Open FITS file GFits file(expression); // Log selected FITS file - if (logExplicit()) { - log << std::endl; - log.header1("FITS file content after selection"); - log << file << std::endl; - } + log_header1(EXPLICIT, "FITS file content after selection"); + log_string(EXPLICIT, file.print(m_chatter)); // Check if we have an events HDU if (!file.contains(evtname)) { @@ -1022,7 +971,8 @@ GEbounds ctselect::set_ebounds(GCTAObservation* obs, const GEbounds& ebounds) co // Get CTA IRF respsonse pointer and throw exception if we don't // find IRF - const GCTAResponseIrf* rsp = dynamic_cast(obs->response()); + const GCTAResponseIrf* rsp = dynamic_cast + (obs->response()); if (rsp == NULL) { std::string msg = "No IRF response attached to given observation."; throw GException::invalid_value(G_SET_EBOUNDS, msg); @@ -1268,10 +1218,7 @@ void ctselect::save_fits(void) } // Log filename - if (logTerse()) { - log << gammalib::parformat("Event list file"); - log << outfile << std::endl; - } + log_value(NORMAL, "Event list file", outfile); // Save event list save_event_list(obs, m_infiles[0], m_evtname[0], m_gtiname[0], @@ -1300,19 +1247,10 @@ void ctselect::save_xml(void) { // Get output filename and prefix m_outobs = (*this)["outobs"].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_outobs.substr(m_outobs.length()-4,4)); - if (suffix != ".xml") { - log << "*** WARNING: Name of observation definition output file \""+ - m_outobs+"\"" << std::endl; - log << "*** WARNING: does not terminate with \".xml\"." << std::endl; - log << "*** WARNING: This is not an error, but might be misleading." - " It is recommended" << std::endl; - log << "*** WARNING: to use the suffix \".xml\" for observation" - " definition files." << std::endl; - } + log_string(TERSE, warn_xml_suffix(m_outobs)); // Loop over all observation in the container for (int i = 0; i < m_obs.size(); ++i) { @@ -1337,10 +1275,7 @@ void ctselect::save_xml(void) outfile += "["+m_evtname[i]+"]"; // Log filename - if (logTerse()) { - log << gammalib::parformat("Event list file"); - log << outfile << std::endl; - } + log_value(NORMAL, "Event list file", outfile); // Store output file name in observation obs->eventfile(outfile); @@ -1351,10 +1286,8 @@ void ctselect::save_xml(void) } // endfor: looped over observations - // Dump filename - if (logTerse()) { - log << "Save \""+m_outobs+"\"" << std::endl; - } + // 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); @@ -1421,7 +1354,8 @@ void ctselect::save_event_list(const GCTAObservation* obs, std::string outevt = evtname; std::string outgti = gtiname; if (outname.has_extname()) { - std::vector extnames = gammalib::split(outname.extname(), ";"); + std::vector extnames = + gammalib::split(outname.extname(), ";"); if (extnames.size() > 0) { std::string extname = gammalib::strip_whitespace(extnames[0]); if (!extname.empty()) { diff --git a/src/ctselect/ctselect.hpp b/src/ctselect/ctselect.hpp index 59c7a48..e73d8c3 100644 --- a/src/ctselect/ctselect.hpp +++ b/src/ctselect/ctselect.hpp @@ -36,7 +36,7 @@ /* __Definitions _________________________________________________________ */ #define CTSELECT_NAME "ctselect" -#define CTSELECT_VERSION "1.0.0" +#define CTSELECT_VERSION "1.2.0" /***********************************************************************//** @@ -102,6 +102,7 @@ protected: double m_emax; //!< Upper energy std::string m_expr; //!< Selection expression std::string m_usethres; //!< Energy threshold type + GChatter m_chatter; //!< Chattiness // Protected members GObservations m_obs; //!< Observations container diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index fb29760..397a3ef 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -221,64 +221,34 @@ void ctskymap::run(void) // Get parameters get_parameters(); - // 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 input observation container into logger + log_observations(NORMAL, m_obs, "Input observation"); - // Write header - if (logTerse()) { - log << std::endl; - if (m_obs.size() > 1) { - log.header1("Map observations"); - } - else { - log.header1("Map observation"); - } - } + // Write header into logger + log_header1(TERSE, gammalib::number("Map observation", m_obs.size())); // Loop over all observation in the container for (int i = 0; i < m_obs.size(); ++i) { - // 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() +")"; - } - log.header3(header); - } + // Write header for the current observation + log_header3(TERSE, get_obs_header(m_obs[i])); // Get CTA observation GCTAObservation* obs = dynamic_cast(m_obs[i]); // Skip observation if it's not CTA if (obs == NULL) { - if (logTerse()) { - log << " Skipping "; - log << m_obs[i]->instrument(); - log << " observation" << std::endl; - } + 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") { - if (logTerse()) { - log << " Skipping binned "; - log << obs->instrument(); - log << " observation" << std::endl; - } + std::string msg = " Skipping binned "+obs->instrument()+ + " observation"; + log_string(NORMAL, msg); continue; } @@ -309,22 +279,16 @@ void ctskymap::run(void) void ctskymap::save(void) { // Write header - if (logTerse()) { - log << std::endl; - log.header1("Save sky map"); - } + log_header1(TERSE, "Save sky map"); // Get sky map filename m_outmap = (*this)["outmap"].filename(); - // Save only if filename is non-empty - if (!m_outmap.is_empty()) { + // Determine whether the map is empty + bool map_is_empty = ((m_skymap.nx() == 0) || (m_skymap.ny() == 0)); - // Log filename - if (logTerse()) { - log << gammalib::parformat("Sky map file"); - log << m_outmap.url() << std::endl; - } + // Save sky map if filename and the map are not empty + if (!m_outmap.is_empty() && !map_is_empty) { // Create empty FITS file GFits fits; @@ -337,6 +301,13 @@ void ctskymap::save(void) } + // Write into logger what has been done + std::string fname = (m_outmap.is_empty()) ? "NONE" : m_outmap.url(); + if (map_is_empty) { + fname.append(" (map is empty, no file created)"); + } + log_value(NORMAL, "Sky map file", fname); + // Return return; } @@ -349,11 +320,8 @@ void ctskymap::save(void) ***************************************************************************/ void ctskymap::publish(const std::string& name) { - // Write header - if (logTerse()) { - log << std::endl; - log.header1("Publish sky map"); - } + // Write header into logger + log_header1(TERSE, "Publish sky map"); // Set default name is user name is empty std::string user_name(name); @@ -361,11 +329,8 @@ void ctskymap::publish(const std::string& name) user_name = CTSKYMAP_NAME; } - // Log filename - if (logTerse()) { - log << gammalib::parformat("Publish sky map"); - log << user_name << std::endl; - } + // Write sky map name into logger + log_value(NORMAL, "Sky map name", user_name); // Publish sky map m_skymap.publish(user_name); @@ -375,6 +340,60 @@ void ctskymap::publish(const std::string& name) } +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctskymap::init_members(void) +{ + // Initialise members + m_obs.clear(); + m_skymap.clear(); + m_emin = 0.0; + m_emax = 0.0; + m_publish = false; + m_chatter = static_cast(2); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctskymap::copy_members(const ctskymap& app) +{ + // Copy attributes + m_obs = app.m_obs; + m_skymap = app.m_skymap; + m_emin = app.m_emin; + m_emax = app.m_emax; + m_publish = app.m_publish; + m_chatter = app.m_chatter; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctskymap::free_members(void) +{ + // Return + return; +} + + /***********************************************************************//** * @brief Get application parameters * @@ -408,6 +427,7 @@ void ctskymap::get_parameters(void) m_emin = (*this)["emin"].real(); m_emax = (*this)["emax"].real(); m_publish = (*this)["publish"].boolean(); + m_chatter = static_cast((*this)["chatter"].integer()); // Read ahead parameters if (read_ahead()) { @@ -436,131 +456,66 @@ void ctskymap::get_parameters(void) ***************************************************************************/ 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); + } - // 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); + + // Initialise binning statistics + int num_outside_map = 0; + int num_outside_erange = 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]; + + // Skip if energy is out of range + if (event->energy() < emin || event->energy() > emax) { + num_outside_erange++; + continue; } - // Setup energy range covered by data - GEnergy emin; - GEnergy emax; - GEbounds ebds; - emin.TeV(m_emin); - emax.TeV(m_emax); - - // Initialise binning statistics - int num_outside_map = 0; - int num_outside_erange = 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]; - - // 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->dir(); - GSkyPixel pixel = m_skymap.dir2pix(dir); - - // Skip if pixel is out of range - if (pixel.x() < -0.5 || pixel.x() > (m_skymap.nx() - 0.5) || - pixel.y() < -0.5 || pixel.y() > (m_skymap.ny() - 0.5)) { - num_outside_map++; - continue; - } - - // Fill event in skymap - m_skymap(pixel, 0) += 1.0; - num_in_map++; - - } // endfor: looped over all events - - // Log binning results - if (logTerse()) { - log << gammalib::parformat("Events in list"); - log << obs->events()->size() << std::endl; - log << gammalib::parformat("Events in map"); - log << num_in_map << std::endl; - log << gammalib::parformat("Events outside map area"); - log << num_outside_map << std::endl; - log << gammalib::parformat("Events outside energies"); - log << num_outside_erange << std::endl; - } + // Determine sky pixel + GCTAInstDir* inst = (GCTAInstDir*)&(event->dir()); + GSkyDir dir = inst->dir(); + GSkyPixel pixel = m_skymap.dir2pix(dir); - // Log map - if (logVerbose()) { - log.indent(1); - log << m_skymap << std::endl; - log.indent(0); + // Skip if pixel is out of range + if (pixel.x() < -0.5 || pixel.x() > (m_skymap.nx() - 0.5) || + pixel.y() < -0.5 || pixel.y() > (m_skymap.ny() - 0.5)) { + num_outside_map++; + continue; } - } // endif: observation was valid + // Fill event in skymap + m_skymap(pixel, 0) += 1.0; + num_in_map++; - // Return - return; -} + } // endfor: looped over all events + // Log binning results + log_value(NORMAL, "Events in list", obs->events()->size()); + log_value(NORMAL, "Events in map", num_in_map); + log_value(NORMAL, "Events outside map area", num_outside_map); + log_value(NORMAL, "Events outside energies", num_outside_erange); -/*========================================================================== - = = - = Private methods = - = = - ==========================================================================*/ + // Write sky map into header + log_header2(EXPLICIT, "Sky map"); + log_string(EXPLICIT, m_skymap.print(m_chatter)); -/***********************************************************************//** - * @brief Initialise class members - ***************************************************************************/ -void ctskymap::init_members(void) -{ - // Initialise members - m_obs.clear(); - m_skymap.clear(); - m_emin = 0.0; - m_emax = 0.0; - m_publish = false; - - // Return - return; -} - - -/***********************************************************************//** - * @brief Copy class members - * - * @param[in] app Application. - ***************************************************************************/ -void ctskymap::copy_members(const ctskymap& app) -{ - // Copy attributes - m_obs = app.m_obs; - m_skymap = app.m_skymap; - m_emin = app.m_emin; - m_emax = app.m_emax; - m_publish = app.m_publish; - - // Return - return; -} - - -/***********************************************************************//** - * @brief Delete class members - ***************************************************************************/ -void ctskymap::free_members(void) -{ // Return return; } diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp index 5f6498b..1636b80 100644 --- a/src/ctskymap/ctskymap.hpp +++ b/src/ctskymap/ctskymap.hpp @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTSKYMAP_NAME "ctskymap" -#define CTSKYMAP_VERSION "1.1.0" +#define CTSKYMAP_VERSION "1.2.0" /***********************************************************************//** @@ -77,6 +77,7 @@ protected: double m_emin; //!< Minimum energy (TeV) double m_emax; //!< Maximum energy (TeV) bool m_publish; //!< Publish sky map? + GChatter m_chatter; //!< Chattiness // Protected members GObservations m_obs; //!< Observation container diff --git a/src/cttsmap/cttsmap.cpp b/src/cttsmap/cttsmap.cpp index 716a8e5..058b2a8 100644 --- a/src/cttsmap/cttsmap.cpp +++ b/src/cttsmap/cttsmap.cpp @@ -215,36 +215,15 @@ void cttsmap::run(void) // Get task parameters get_parameters(); - // Set energy dispersion flag for all CTA observations and save old + // Set energy dispersion flags of all CTA observations and save old // values in save_edisp vector - std::vector save_edisp; - save_edisp.assign(m_obs.size(), false); - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - save_edisp[i] = obs->response()->apply_edisp(); - obs->response()->apply_edisp(m_apply_edisp); - } - } + std::vector save_edisp = set_edisp(m_obs, m_apply_edisp); - // 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 input observation container into logger + log_observations(NORMAL, m_obs, "Input observation"); // Write header - if (logTerse()) { - log << std::endl; - log.header1("Initialise TS map"); - } - + log_header1(TERSE, "Initialise TS map"); // Get bins to be computed // To split the computation on several jobs, @@ -270,10 +249,7 @@ void cttsmap::run(void) if (m_logL0 == 0.0) { // Write header - if (logTerse()) { - log << std::endl; - log.header1("Compute NULL Hypothesis for TS computation"); - } + log_header1(TERSE, "Compute NULL Hypothesis for TS computation"); // Compute likelihood without the test source m_obs.models(models); @@ -283,10 +259,7 @@ void cttsmap::run(void) } // Write header - if (logTerse()) { - log << std::endl; - log.header1("Generate TS map"); - } + log_header1(TERSE, "Generate TS map"); // Loop over grid positions for (int i = binmin; i < binmax; ++i) { @@ -295,12 +268,8 @@ void cttsmap::run(void) GSkyDir bincentre = m_tsmap.inx2dir(i); // Header for verbose logging - if (logExplicit()) { - std::string msg = "Computing TS for bin number "+gammalib::str(i)+ - " at "+bincentre.print(); - log << std::endl; - log.header2(msg); - } + log_header2(EXPLICIT, "Computing TS for bin number "+gammalib::str(i)+ + " at "+bincentre.print()); // Add test source at current bin position if (m_testsource != NULL) { @@ -360,13 +329,8 @@ void cttsmap::run(void) // Bring models to initial state m_obs.models(models_orig); - // Restore energy dispersion flag for all CTA observations - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - obs->response()->apply_edisp(save_edisp[i]); - } - } + // Restore energy dispersion flags of all CTA observations + restore_edisp(m_obs, save_edisp); // Optionally publish TS map if (m_publish) { @@ -386,22 +350,16 @@ void cttsmap::run(void) void cttsmap::save(void) { // Write header - if (logTerse()) { - log << std::endl; - log.header1("Save TS map"); - } + log_header1(TERSE, "Save TS map"); // Get TS map filename m_outmap = (*this)["outmap"].filename(); - // Save only if filename is non-empty - if (!m_outmap.is_empty()) { + // Determine whether the TS map is empty + bool map_is_empty = ((m_tsmap.nx() == 0) || (m_tsmap.ny() == 0)); - // Log filename - if (logTerse()) { - log << gammalib::parformat("TS map file"); - log << m_outmap.url() << std::endl; - } + // Save only if filename is non-empty + if (!m_outmap.is_empty() && !map_is_empty) { // Create fits file GFits fits; @@ -428,6 +386,13 @@ void cttsmap::save(void) } // endif: TS map filename was valid + // Write into logger what has been done + std::string fname = (m_outmap.is_empty()) ? "NONE" : m_outmap.url(); + if (map_is_empty) { + fname.append(" (TS map is empty, no file created)"); + } + log_value(NORMAL, "TS map file", fname); + // Return return; } @@ -440,11 +405,8 @@ void cttsmap::save(void) ***************************************************************************/ void cttsmap::publish(const std::string& name) { - // Write header - if (logTerse()) { - log << std::endl; - log.header1("Publish TS map"); - } + // Write header into logger + log_header1(TERSE, "Publish TS map"); // Set default name is user name is empty std::string user_name(name); @@ -452,11 +414,8 @@ void cttsmap::publish(const std::string& name) user_name = CTTSMAP_NAME; } - // Log filename - if (logTerse()) { - log << gammalib::parformat("Publish TS map"); - log << user_name << std::endl; - } + // Write TS map name into logger + log_value(NORMAL, "TS map name", user_name); // Publish TS map m_tsmap.publish(user_name); diff --git a/src/cttsmap/cttsmap.hpp b/src/cttsmap/cttsmap.hpp index 49ff6cf..9a21805 100644 --- a/src/cttsmap/cttsmap.hpp +++ b/src/cttsmap/cttsmap.hpp @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTTSMAP_NAME "cttsmap" -#define CTTSMAP_VERSION "1.1.0" +#define CTTSMAP_VERSION "1.2.0" /***********************************************************************//** diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index c66c6f1..df5fcea 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -1543,6 +1543,49 @@ GEnergies ctool::insert_energy_boundaries(const GEnergies& energies, /***********************************************************************//** + * @brief Determine the counts cube layer usage. + * + * @param[in] ebounds_cube Energy boundaries of the counts cube. + * @param[in] ebounds_list Energy boundaries of the event list. + * @return Vector of usage flags. + * + * Determines a vector of counts cube layer usage flags that signal whether + * an event should actually be filled in the counts cube or not. This makes + * sure that no partially filled bins will exist in the counts cube. + ***************************************************************************/ +std::vector ctool::cube_layer_usage(const GEbounds& cube_ebounds, + const GEbounds& list_ebounds) const +{ + // Set energy margin + const GEnergy energy_margin(0.01, "GeV"); + + // Initialise usage vector + std::vector usage(cube_ebounds.size(), true); + + // Loop over all energy bins of the cube + for (int iebin = 0; iebin < cube_ebounds.size(); ++iebin) { + + // If the counts cube energy bin is fully contained within the + // energy boundaries of the event list the signal usage of this + // counts cube bin. Partially overlapping energy bins are signaled + // for non-usage, which avoids having partially filled bins in the + // counts cube. Some margin is applied that effectively reduces the + // width of the counts cube energy bin. This should cope with any + // rounding errors that come from reading and writing the energy + // boundary information to a FITS file. + usage[iebin] = list_ebounds.contains(cube_ebounds.emin(iebin) + + energy_margin, + cube_ebounds.emax(iebin) - + energy_margin); + + } // endfor: looped over all energy bins + + // Return usage + return usage; +} + + +/***********************************************************************//** * @brief Set warning string if there are too few energies * * @param[in] energies Energies. @@ -1584,6 +1627,37 @@ std::string ctool::warn_too_few_energies(const GEnergies& energies) const } +/***********************************************************************//** + * @brief Set warning string if file has no .xml suffix + * + * @param[in] filename Filename. + * @return Warning string. + * + * Sets a warning string if the filename has not ".xml" suffix. + ***************************************************************************/ +std::string ctool::warn_xml_suffix(const GFilename& filename) const +{ + // Initialise warning + std::string warning; + + // Get filename suffix + std::string fname = std::string(filename); + std::string suffix = gammalib::tolower(fname.substr(fname.length()-4,4)); + + // Warn if there are not enough energy bins + if (suffix != ".xml") { + warning.append("\nWARNING: Name of observation definition output file " + "\""+filename+"\" does not terminate with \".xml\"." + "\n This is not an error, but may be misleading. " + "It is recommended to use the suffix" + "\n \".xml\" for observation definition files."); + } + + // Return warning + return warning; +} + + /*========================================================================== = = = Private methods for SWIG = diff --git a/src/support/ctool.hpp b/src/support/ctool.hpp index 7178f0a..427a24e 100644 --- a/src/support/ctool.hpp +++ b/src/support/ctool.hpp @@ -105,9 +105,12 @@ protected: std::string get_obs_header(const GObservation* obs); GEnergies insert_energy_boundaries(const GEnergies& energies, const GCTAObservation& obs); + std::vector cube_layer_usage(const GEbounds& cube_ebounds, + const GEbounds& list_ebounds) const; // Protected warning strings std::string warn_too_few_energies(const GEnergies& energies) const; + std::string warn_xml_suffix(const GFilename& filename) const; // Protected members bool m_read_ahead; //!< Read ahead output parameters diff --git a/test/test_ctmapcube.py b/test/test_ctmapcube.py index e5604db..cc45c2c 100644 --- a/test/test_ctmapcube.py +++ b/test/test_ctmapcube.py @@ -18,6 +18,7 @@ # along with this program. If not, see . # # ========================================================================== +import os import gammalib import ctools from testing import test @@ -85,7 +86,7 @@ class Test(test): 'Check successful execution from command line') # Check map cube - self._check_cube('ctmapcube_cmd1.fits') + self._check_result_file('ctmapcube_cmd1.fits') # Setup ctmapcube command cmd = ctmapcube+' inmodel="events_that_do_not_exist.fits"'+ \ @@ -107,10 +108,26 @@ class Test(test): """ Test ctmapcube from Python """ - # Set-up ctmapcube + # Allocate ctmapcube mapcube = ctools.ctmapcube() + + # Check that empty ctmapcube tool holds an map cube that has no + # energy bins + self._check_cube(mapcube.mapcube(), nmaps=0, npixels=0) + + # Check that saving does not nothing + mapcube['outcube'] = 'ctmapcube_py0.fits' + mapcube['logfile'] = 'ctmapcube_py0.log' + mapcube.logFileOpen() + mapcube.save() + self.test_assert(not os.path.isfile('ctmapcube_py0.fits'), + 'Check that no map cube has been created') + + # Check that clearing does not lead to an exception or segfault + mapcube.clear() + + # Now set ctmapcube parameters mapcube['inmodel'] = self._model - mapcube['outcube'] = 'ctmapcube_py1.fits' mapcube['ebinalg'] = 'LOG' mapcube['emin'] = 0.1 mapcube['emax'] = 100.0 @@ -122,6 +139,7 @@ class Test(test): mapcube['proj'] = 'CAR' mapcube['xref'] = 83.63 mapcube['yref'] = 22.01 + mapcube['outcube'] = 'ctmapcube_py1.fits' mapcube['logfile'] = 'ctmapcube_py1.log' mapcube['chatter'] = 2 @@ -133,31 +151,76 @@ class Test(test): mapcube.save() # Check map cube - self._check_cube('ctmapcube_py1.fits') + self._check_result_file('ctmapcube_py1.fits') + + # Copy ctmapcube tool + cpy_mapcube = mapcube.copy() + + # Check map cube of ctmapcube copy + self._check_cube(cpy_mapcube.mapcube()) + + # Execute copy of ctmapcube tool again, now with a higher chatter + # level than before + cpy_mapcube['outcube'] = 'ctmapcube_py2.fits' + cpy_mapcube['logfile'] = 'ctmapcube_py2.log' + cpy_mapcube['publish'] = True + cpy_mapcube['chatter'] = 3 + cpy_mapcube.logFileOpen() # Needed to get a new log file + cpy_mapcube.execute() + + # Check result file + self._check_result_file('ctmapcube_py2.fits') + + # Now clear copy of ctmapcube tool + cpy_mapcube.clear() + + # Check that the cleared copy has also cleared the map cube + self._check_cube(cpy_mapcube.mapcube(), nmaps=0, npixels=0) # Return return - # Check map cube - def _check_cube(self, filename): + # Check result file + def _check_result_file(self, filename): """ Check content of map cube Parameters ---------- filename : str - Map cube file name. + Map cube file name """ # Load map cube cube = gammalib.GModelSpatialDiffuseCube(filename) - # Test map cube - #self.test_value(cube.maps(), 10, "10 maps") - #self.test_value(cube.pixels(), 10000, "10000 map pixels") - # The map cube is not loaded by default !!!! We should add the method so that - # the attributes are known !!!! - self.test_value(cube.maps(), 0, '10 maps') - self.test_value(cube.pixels(), 0, '10000 map pixels') + # Check map cube + self._check_cube(cube) + + # Return + return + + # Check map cube + def _check_cube(self, cube, nmaps=10, npixels=10000): + """ + Check map cube + + Parameters + ---------- + cube : `~gammalib.GModelSpatialDiffuseCube` + Map cube + nmaps : int, optional + Number of maps + npixels : int, optional + Number of pixels + """ + # Get energies (this forces loading in case the map cube is not + # loaded) + energies = cube.energies() + self.test_value(len(cube.energies()), nmaps, 'Check number of energies') + + # Check dimensions + self.test_value(cube.maps(), nmaps, 'Check number of maps') + self.test_value(cube.pixels(), npixels, 'Check number of pixels') # Return return diff --git a/test/test_ctmodel.py b/test/test_ctmodel.py index 7c28c60..6534525 100755 --- a/test/test_ctmodel.py +++ b/test/test_ctmodel.py @@ -18,6 +18,7 @@ # along with this program. If not, see . # # ========================================================================== +import os import gammalib import ctools from testing import test @@ -71,8 +72,8 @@ class Test(test): ' inmodel="'+self._model+'"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ ' rad=5.0 ra=83.63 dec=22.01 tmin=0.0 tmax=1800'+ \ - ' emin=0.1 emax=100.0 enumbins=20 nxpix=200 nypix=200'+ \ - ' binsz=0.02 coordsys="CEL" proj="CAR"'+ \ + ' emin=0.1 emax=100.0 enumbins=10 nxpix=40 nypix=40'+ \ + ' binsz=0.1 coordsys="CEL" proj="CAR"'+ \ ' xref=83.63 yref=22.01'+ \ ' logfile="ctmodel_cmd1.log" chatter=1' @@ -94,8 +95,8 @@ class Test(test): ' inmodel="model_that_does_not_exist.xml"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ ' rad=5.0 ra=83.63 dec=22.01 tmin=0.0 tmax=1800'+ \ - ' emin=0.1 emax=100.0 enumbins=20 nxpix=200 nypix=200'+ \ - ' binsz=0.02 coordsys="CEL" proj="CAR"'+ \ + ' emin=0.1 emax=100.0 enumbins=10 nxpix=40 nypix=40'+ \ + ' binsz=0.1 coordsys="CEL" proj="CAR"'+ \ ' xref=83.63 yref=22.01'+ \ ' logfile="ctmodel_cmd2.log" chatter=1' @@ -111,10 +112,25 @@ class Test(test): """ Test ctmodel from Python """ - # Set-up ctmodel from scratch + # Allocate ctmodel model = ctools.ctmodel() + + # Check that empty ctmodel tool holds an empty model cube + self._check_cube(model.cube(), nx=0, ny=0, nebins=0) + + # Check that saving does not nothing + model['outcube'] = 'ctmodel_py0.fits' + model['logfile'] = 'ctmodel_py0.log' + model.logFileOpen() + model.save() + self.test_assert(not os.path.isfile('ctmodel_py0.fits'), + 'Check that no model cube has been created') + + # Check that clearing does not lead to an exception or segfault + model.clear() + + # Now set ctmodel parameters model['incube'] = 'NONE' - model['outcube'] = 'ctmodel_py1.fits' model['inmodel'] = self._model model['inobs'] = 'NONE' model['expcube'] = 'NONE' @@ -128,15 +144,16 @@ class Test(test): model['tmin'] = 0 model['tmax'] = 1800 model['emin'] = 0.1 - model['emax'] = 100 - model['enumbins'] = 20 - model['nxpix'] = 200 - model['nypix'] = 200 - model['binsz'] = 0.02 + model['emax'] = 100.0 + model['enumbins'] = 10 + model['nxpix'] = 40 + model['nypix'] = 40 + model['binsz'] = 0.1 model['coordsys'] = 'CEL' model['proj'] = 'CAR' model['xref'] = 83.63 model['yref'] = 22.01 + model['outcube'] = 'ctmodel_py1.fits' model['logfile'] = 'ctmodel_py1.log' model['chatter'] = 2 @@ -148,6 +165,30 @@ class Test(test): # Check result file self._check_result_file('ctmodel_py1.fits') + # Copy ctmodel tool + cpy_model = model.copy() + + # Check model cube of ctmodel copy + self._check_cube(cpy_model.cube()) + + # Execute copy of ctmodel tool again, now with a higher chatter + # level than before + cpy_model['outcube'] = 'ctmodel_py2.fits' + cpy_model['logfile'] = 'ctmodel_py2.log' + cpy_model['publish'] = True + cpy_model['chatter'] = 3 + cpy_model.logFileOpen() # Needed to get a new log file + cpy_model.execute() + + # Check result file + self._check_result_file('ctmodel_py2.fits') + + # Now clear copy of ctmodel tool + cpy_model.clear() + + # Check that the cleared copy has also cleared the model cube + self._check_cube(cpy_model.cube(), nx=0, ny=0, nebins=0) + # Return return @@ -155,20 +196,41 @@ class Test(test): def _check_result_file(self, filename): """ Check result file + + Parameters + ---------- + filename : str + Model cube file name """ # Open result file - fits = gammalib.GFits(filename) - - # Get HDUs - cube = fits['Primary'] - ebounds = fits['EBOUNDS'] - gti = fits['GTI'] - - # Check dimensions - self.test_value(cube.naxis(), 3, 'Check for 3 cube dimensions') - self.test_value(cube.naxes(0), 200, 'Check for 200 pixels in X') - self.test_value(cube.naxes(1), 200, 'Check for 200 pixels in Y') - self.test_value(cube.naxes(2), 20, 'Check for 20 pixels in Z') + cube = gammalib.GCTAEventCube(filename) + + # Check cube + self._check_cube(cube) + + # Return + return + + # Check model cube + def _check_cube(self, cube, nx=40, ny=40, nebins=10): + """ + Check model cube + + Parameters + ---------- + cube : `~gammalib.GCTAEventCube` + Model cube + nx : int, optional + Number of pixels in X direction + ny : int, optional + Number of pixels in Y direction + nebins : int, optional + Number of energy bins + """ + # Check number of bins per axis + self.test_value(cube.nx(), nx, 'Check number of X pixels') + self.test_value(cube.ny(), ny, 'Check number of Y pixels') + self.test_value(cube.ebins(), nebins, 'Check number of energy bins') # Return return