Iftach Sadeh / ctools

Commit 0622a2cb4b42dce0642abb4cc4d5eee7e4a8e100

Authored by Jürgen Knödlseder
1 parent 115eb8a9

Use new GApplication logger methods and enhance ctool unit tests (#1828)

Inline Side-by-side
Showing 19 changed files with 690 additions and 780 deletions
pyext/ctmapcube.i
  View file @0622a2c
... ... @@ -44,11 +44,12 @@ public:
44 44 virtual ~ctmapcube(void);
45 45  
46 46 // Methods
47   - void clear(void);
48   - void run(void);
49   - void save(void);
50   - void publish(const std::string& name = "");
51   - void models(const GModels& models);
  47 + void clear(void);
  48 + void run(void);
  49 + void save(void);
  50 + void publish(const std::string& name = "");
  51 + const GModelSpatialDiffuseCube& mapcube(void) const;
  52 + void models(const GModels& models);
52 53 };
53 54  
54 55  
... ...
pyext/ctool.i
  View file @0622a2c
... ... @@ -84,26 +84,31 @@ public:
84 84 virtual void execute(void);
85 85  
86 86 // Make methods private in Python by prepending an underscore
87   - %rename(_read_ahead) read_ahead;
88   - %rename(_time_reference) time_reference;
89   - %rename(_get_observations) get_observations;
90   - %rename(_setup_observations) setup_observations;
91   - %rename(_create_ebounds) create_ebounds;
92   - %rename(_create_map) create_map;
93   - %rename(_create_cube) create_cube;
94   - %rename(_create_cta_obs) create_cta_obs;
95   - %rename(_require_inobs) require_inobs;
96   - %rename(_require_inobs_nolist) require_inobs_nolist;
97   - %rename(_require_inobs_nocube) require_inobs_nocube;
98   - %rename(_log_observations) log_observations;
99   - %rename(_log_models) log_models;
100   - %rename(_set_response) set_response;
101   - %rename(_set_obs_response) set_obs_response;
102   - %rename(_set_obs_bounds) set_obs_bounds;
103   - %rename(_get_mean_pointing) get_mean_pointing;
104   - %rename(_get_current_rss) get_current_rss;
105   - %rename(_get_obs_header) get_obs_header;
106   - %rename(_warn_too_few_energies) warn_too_few_energies;
  87 + %rename(_read_ahead) read_ahead;
  88 + %rename(_time_reference) time_reference;
  89 + %rename(_get_observations) get_observations;
  90 + %rename(_setup_observations) setup_observations;
  91 + %rename(_create_ebounds) create_ebounds;
  92 + %rename(_create_map) create_map;
  93 + %rename(_create_cube) create_cube;
  94 + %rename(_create_cta_obs) create_cta_obs;
  95 + %rename(_require_inobs) require_inobs;
  96 + %rename(_require_inobs_nolist) require_inobs_nolist;
  97 + %rename(_require_inobs_nocube) require_inobs_nocube;
  98 + %rename(_log_observations) log_observations;
  99 + %rename(_log_models) log_models;
  100 + %rename(_set_response) set_response;
  101 + %rename(_set_edisp) set_edisp;
  102 + %rename(_restore_edisp) restore_edisp;
  103 + %rename(_set_obs_response) set_obs_response;
  104 + %rename(_set_obs_bounds) set_obs_bounds;
  105 + %rename(_get_mean_pointing) get_mean_pointing;
  106 + %rename(_get_current_rss) get_current_rss;
  107 + %rename(_get_obs_header) get_obs_header;
  108 + %rename(_insert_energy_boundaries) insert_energy_boundaries;
  109 + %rename(_cube_layer_usage) cube_layer_usage;
  110 + %rename(_warn_too_few_energies) warn_too_few_energies;
  111 + %rename(_warn_xml_suffix) warn_xml_suffix;
107 112  
108 113 // Protected methods
109 114 const bool& read_ahead(void) const;
... ... @@ -131,16 +136,23 @@ public:
131 136 const std::string& what = "Model");
132 137  
133 138 // Protected support methods
134   - void set_response(GObservations& obs);
135   - void set_obs_response(GCTAObservation* obs);
136   - void set_obs_bounds(GObservations& obs);
137   - GSkyDir get_mean_pointing(const GObservations& obs);
138   - size_t get_current_rss(void);
139   - std::string get_obs_header(const GObservation* obs);
140   - GEnergies insert_energy_boundaries(const GEnergies& energies,
141   - const GCTAObservation& obs);
  139 + void set_response(GObservations& obs);
  140 + std::vector<bool> set_edisp(GObservations& obs, const bool& edisp) const;
  141 + void restore_edisp(GObservations& obs,
  142 + const std::vector<bool>& edisp) const;
  143 + void set_obs_response(GCTAObservation* obs);
  144 + void set_obs_bounds(GObservations& obs);
  145 + GSkyDir get_mean_pointing(const GObservations& obs);
  146 + size_t get_current_rss(void);
  147 + std::string get_obs_header(const GObservation* obs);
  148 + GEnergies insert_energy_boundaries(const GEnergies& energies,
  149 + const GCTAObservation& obs);
  150 + std::vector<bool> cube_layer_usage(const GEbounds& cube_ebounds,
  151 + const GEbounds& list_ebounds) const;
  152 +
142 153 // Protected warning strings
143 154 std::string warn_too_few_energies(const GEnergies& energies) const;
  155 + std::string warn_xml_suffix(const GFilename& filename) const;
144 156 };
145 157  
146 158  
... ...
src/ctbin/ctbin.cpp
  View file @0622a2c
... ... @@ -254,7 +254,7 @@ void ctbin::run(void)
254 254  
255 255 // Skip observation if we have a binned observation
256 256 if (obs->eventtype() == "CountsCube") {
257   - std::string msg = " Skipping binned "+m_obs[i]->instrument()+
  257 + std::string msg = " Skipping binned "+obs->instrument()+
258 258 " observation";
259 259 log_string(NORMAL, msg);
260 260 continue;
... ... @@ -518,7 +518,7 @@ void ctbin::fill_cube(GCTAObservation* obs)
518 518 }
519 519  
520 520 // Get counts cube usage flags
521   - std::vector<bool> usage = cube_layer_usage(events->ebounds());
  521 + std::vector<bool> usage = cube_layer_usage(m_ebounds, events->ebounds());
522 522  
523 523 // Initialise binning statistics
524 524 int num_outside_roi = 0;
... ... @@ -624,7 +624,7 @@ void ctbin::set_weights(GCTAObservation* obs)
624 624 }
625 625  
626 626 // Get counts cube usage flags
627   - std::vector<bool> usage = cube_layer_usage(events->ebounds());
  627 + std::vector<bool> usage = cube_layer_usage(m_ebounds, events->ebounds());
628 628  
629 629 // Loop over all pixels in counts cube
630 630 for (int pixel = 0; pixel < m_counts.npix(); ++pixel) {
... ... @@ -778,42 +778,3 @@ void ctbin::obs_cube(void)
778 778 // Return
779 779 return;
780 780 }
781   -
782   -
783   -/***********************************************************************//**
784   - * @brief Determine the counts cube layer usage.
785   - *
786   - * @param[in] ebounds Energy boundaries of the event list.
787   - * @return Vector of usage flags.
788   - *
789   - * Determines a vector of counts cube layer usage flags that signal whether
790   - * an event should actually be filled in the counts cube or not. This makes
791   - * sure that no partially filled bins will exist in the counts cube.
792   - ***************************************************************************/
793   -std::vector<bool> ctbin::cube_layer_usage(const GEbounds& ebounds) const
794   -{
795   - // Set energy margin
796   - const GEnergy energy_margin(0.01, "GeV");
797   -
798   - // Initialise usage vector
799   - std::vector<bool> usage(m_ebounds.size(), true);
800   -
801   - // Loop over all energy bins of the cube
802   - for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
803   -
804   - // If the counts cube energy bin is fully contained within the
805   - // energy boundaries of the event list the signal usage of this
806   - // counts cube bin. Partially overlapping energy bins are signaled
807   - // for non-usage, which avoids having partially filled bins in the
808   - // counts cube. Some margin is applied that effectively reduces the
809   - // width of the counts cube energy bin. This should cope with any
810   - // rounding errors that come from reading and writing the energy
811   - // boundary information to a FITS file.
812   - usage[iebin] = ebounds.contains(m_ebounds.emin(iebin) + energy_margin,
813   - m_ebounds.emax(iebin) - energy_margin);
814   -
815   - } // endfor: looped over all energy bins
816   -
817   - // Return usage
818   - return usage;
819   -}
... ...
src/ctbin/ctbin.hpp
  View file @0622a2c
... ... @@ -85,14 +85,13 @@ public:
85 85  
86 86 protected:
87 87 // Protected methods
88   - void init_members(void);
89   - void copy_members(const ctbin& app);
90   - void free_members(void);
91   - void get_parameters(void);
92   - void fill_cube(GCTAObservation* obs);
93   - void set_weights(GCTAObservation* obs);
94   - void obs_cube(void);
95   - std::vector<bool> cube_layer_usage(const GEbounds& ebounds) const;
  88 + void init_members(void);
  89 + void copy_members(const ctbin& app);
  90 + void free_members(void);
  91 + void get_parameters(void);
  92 + void fill_cube(GCTAObservation* obs);
  93 + void set_weights(GCTAObservation* obs);
  94 + void obs_cube(void);
96 95  
97 96 // User parameters
98 97 GFilename m_outcube; //!< Output counts map file name
... ...
src/ctmapcube/ctmapcube.cpp
  View file @0622a2c
... ... @@ -192,16 +192,16 @@ void ctmapcube::run(void)
192 192 // Get task parameters
193 193 get_parameters();
194 194  
195   - // Write models into logger
196   - if (logTerse()) {
197   - log << std::endl;
198   - log.header1(gammalib::number("Model", m_models.size()));
199   - log << m_models << std::endl;
200   - }
  195 + // Write input model container into logger
  196 + log_models(NORMAL, m_models, "Input model");
201 197  
202 198 // Generate map cube
203 199 create_cube();
204 200  
  201 + // Write result cube into logger
  202 + log_header1(NORMAL, "Map cube");
  203 + log_string(NORMAL, m_cube.print(m_chatter));
  204 +
205 205 // Optionally publish map cube
206 206 if (m_publish) {
207 207 publish();
... ... @@ -220,27 +220,27 @@ void ctmapcube::run(void)
220 220 void ctmapcube::save(void)
221 221 {
222 222 // Write header
223   - if (logTerse()) {
224   - log << std::endl;
225   - log.header1("Save map cube");
226   - }
  223 + log_header1(TERSE, "Save map cube");
227 224  
228 225 // Get map cube filename
229 226 m_outcube = (*this)["outcube"].filename();
230 227  
231   - // Save only if filename is non-empty
232   - if (!m_outcube.is_empty()) {
233   -
234   - // Log filename
235   - if (logTerse()) {
236   - log << gammalib::parformat("Map cube file");
237   - log << m_outcube.url() << std::endl;
238   - }
  228 + // Determine whether the cube is empty
  229 + bool cube_is_empty = ((m_cube.cube().nx() == 0) ||
  230 + (m_cube.cube().ny() == 0) ||
  231 + (m_cube.cube().nmaps() == 0));
239 232  
240   - // Save map cube into FITS file
  233 + // Save map cube if filename and the map cube are not empty
  234 + if (!m_outcube.is_empty() && !cube_is_empty) {
241 235 m_cube.save(m_outcube, clobber());
  236 + }
242 237  
  238 + // Write into logger what has been done
  239 + std::string fname = (m_outcube.is_empty()) ? "NONE" : m_outcube.url();
  240 + if (cube_is_empty) {
  241 + fname.append(" (cube is empty, no file created)");
243 242 }
  243 + log_value(NORMAL, "Map cube file", fname);
244 244  
245 245 // Return
246 246 return;
... ... @@ -256,11 +256,8 @@ void ctmapcube::save(void)
256 256 ***************************************************************************/
257 257 void ctmapcube::publish(const std::string& name)
258 258 {
259   - // Write header
260   - if (logTerse()) {
261   - log << std::endl;
262   - log.header1("Publish map cube");
263   - }
  259 + // Write header into logger
  260 + log_header1(TERSE, "Publish map cube");
264 261  
265 262 // Set default name is user name is empty
266 263 std::string user_name(name);
... ... @@ -268,11 +265,8 @@ void ctmapcube::publish(const std::string&amp; name)
268 265 user_name = CTMAPCUBE_NAME;
269 266 }
270 267  
271   - // Log filename
272   - if (logTerse()) {
273   - log << gammalib::parformat("Map cube");
274   - log << user_name << std::endl;
275   - }
  268 + // Write map cube name into logger
  269 + log_value(NORMAL, "Map cube name", user_name);
276 270  
277 271 // Publish map cube
278 272 m_cube.cube().publish(user_name);
... ... @@ -297,6 +291,7 @@ void ctmapcube::init_members(void)
297 291 m_outcube.clear();
298 292 m_ptsrcsig = 1.0;
299 293 m_publish = false;
  294 + m_chatter = static_cast<GChatter>(2);
300 295  
301 296 // Initialise protected members
302 297 m_models.clear();
... ... @@ -318,6 +313,7 @@ void ctmapcube::copy_members(const ctmapcube&amp; app)
318 313 m_outcube = app.m_outcube;
319 314 m_ptsrcsig = app.m_ptsrcsig;
320 315 m_publish = app.m_publish;
  316 + m_chatter = app.m_chatter;
321 317  
322 318 // Copy protected members
323 319 m_models = app.m_models;
... ... @@ -372,6 +368,7 @@ void ctmapcube::get_parameters(void)
372 368 // Get remaining hidden parameters
373 369 m_ptsrcsig = (*this)["ptsrcsig"].real();
374 370 m_publish = (*this)["publish"].boolean();
  371 + m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
375 372  
376 373 // Read optionally output cube filenames
377 374 if (read_ahead()) {
... ... @@ -462,11 +459,8 @@ GEnergies ctmapcube::create_energies(void)
462 459 ***************************************************************************/
463 460 void ctmapcube::create_cube(void)
464 461 {
465   - // Write header
466   - if (logTerse()) {
467   - log << std::endl;
468   - log.header1("Generate map cube");
469   - }
  462 + // Write header into logger
  463 + log_header1(TERSE, "Generate map cube");
470 464  
471 465 // Loop over all models
472 466 for (int i = 0; i < m_models.size(); ++i) {
... ... @@ -478,10 +472,7 @@ void ctmapcube::create_cube(void)
478 472 if (model == NULL) {
479 473  
480 474 // Signal that model is skipped
481   - if (logNormal()) {
482   - log << gammalib::parformat("Skip model");
483   - log << m_models[i]->name() << std::endl;
484   - }
  475 + log_value(NORMAL, "Skip model", m_models[i]->name());
485 476  
486 477 // Continue
487 478 continue;
... ... @@ -489,10 +480,7 @@ void ctmapcube::create_cube(void)
489 480 }
490 481  
491 482 // Signal that model is used
492   - if (logNormal()) {
493   - log << gammalib::parformat("Use model");
494   - log << m_models[i]->name() << std::endl;
495   - }
  483 + log_value(NORMAL, "Use model", m_models[i]->name());
496 484  
497 485 // Add sky model
498 486 if (dynamic_cast<GModelSpatialPointSource*>(model->spatial()) != NULL) {
... ... @@ -504,13 +492,6 @@ void ctmapcube::create_cube(void)
504 492  
505 493 } // endfor: looped over energies
506 494  
507   - // Log cube
508   - if (logTerse()) {
509   - log << std::endl;
510   - log.header1("Map cube");
511   - log << m_cube << std::endl;
512   - }
513   -
514 495 // Return
515 496 return;
516 497 }
... ...
src/ctmapcube/ctmapcube.hpp
  View file @0622a2c
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTMAPCUBE_NAME "ctmapcube"
37   -#define CTMAPCUBE_VERSION "1.0.0"
  37 +#define CTMAPCUBE_VERSION "1.2.0"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -55,11 +55,12 @@ public:
55 55 ctmapcube& operator=(const ctmapcube& app);
56 56  
57 57 // Methods
58   - void clear(void);
59   - void run(void);
60   - void save(void);
61   - void publish(const std::string& name = "");
62   - void models(const GModels& models);
  58 + void clear(void);
  59 + void run(void);
  60 + void save(void);
  61 + void publish(const std::string& name = "");
  62 + const GModelSpatialDiffuseCube& mapcube(void) const;
  63 + void models(const GModels& models);
63 64  
64 65 protected:
65 66 // Protected methods
... ... @@ -79,6 +80,7 @@ protected:
79 80 GFilename m_outcube; //!< Output map cube filename
80 81 double m_ptsrcsig; //!< Point source sigma (arcmin)
81 82 bool m_publish; //!< Publish map cube?
  83 + GChatter m_chatter; //!< Chattiness
82 84  
83 85 // Protected members
84 86 GModels m_models; //!< Model container
... ... @@ -87,6 +89,18 @@ protected:
87 89  
88 90  
89 91 /***********************************************************************//**
  92 + * @brief Return map cube
  93 + *
  94 + * @return Map cube
  95 + ***************************************************************************/
  96 +inline
  97 +const GModelSpatialDiffuseCube& ctmapcube::mapcube(void) const
  98 +{
  99 + return (m_cube);
  100 +}
  101 +
  102 +
  103 +/***********************************************************************//**
90 104 * @brief Set models
91 105 *
92 106 * @param[in] models Model container.
... ...
src/ctmodel/ctmodel.cpp
  View file @0622a2c
... ... @@ -34,6 +34,7 @@
34 34  
35 35 /* __ Method name definitions ____________________________________________ */
36 36 #define G_GET_PARAMETERS "ctmodel::get_parameters()"
  37 +#define G_FILL_CUBE "ctmodel::fill_cube(GCTAObservation*)"
37 38  
38 39 /* __ Debug definitions __________________________________________________ */
39 40  
... ... @@ -217,68 +218,33 @@ void ctmodel::run(void)
217 218 // Get task parameters
218 219 get_parameters();
219 220  
220   - // Set energy dispersion flag for all CTA observations and save old
  221 + // Set energy dispersion flags of all CTA observations and save old
221 222 // values in save_edisp vector
222   - std::vector<bool> save_edisp;
223   - save_edisp.assign(m_obs.size(), false);
224   - for (int i = 0; i < m_obs.size(); ++i) {
225   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
226   - if (obs != NULL) {
227   - save_edisp[i] = obs->response()->apply_edisp();
228   - obs->response()->apply_edisp(m_apply_edisp);
229   - }
230   - }
  223 + std::vector<bool> save_edisp = set_edisp(m_obs, m_apply_edisp);
231 224  
232   - // Write observation(s) into logger
233   - if (logTerse()) {
234   - log << std::endl;
235   - if (m_obs.size() > 1) {
236   - log.header1("Observations");
237   - }
238   - else {
239   - log.header1("Observation");
240   - }
241   - log << m_obs << std::endl;
242   - }
  225 + // Write input observation container into logger
  226 + log_observations(NORMAL, m_obs, "Input observation");
243 227  
244   - // Write models into logger
245   - if (logTerse()) {
246   - log << std::endl;
247   - log.header1("Models");
248   - log << m_obs.models() << std::endl;
249   - }
  228 + // Write input model container into logger
  229 + log_models(NORMAL, m_obs.models(), "Input model");
250 230  
251 231 // Write header
252   - if (logTerse()) {
253   - log << std::endl;
254   - log.header1("Generate model cube");
255   - }
  232 + log_header1(TERSE, "Generate model cube");
256 233  
257 234 // Loop over all observations in the container
258 235 for (int i = 0; i < m_obs.size(); ++i) {
259 236  
260   - // Write header for observation
261   - if (logTerse()) {
262   - std::string header = m_obs[i]->instrument() + " observation";
263   - if (m_obs[i]->name().length() > 1) {
264   - header += " \"" + m_obs[i]->name() + "\"";
265   - }
266   - if (m_obs[i]->id().length() > 1) {
267   - header += " (id=" + m_obs[i]->id() +")";
268   - }
269   - log.header3(header);
270   - }
  237 + // Write header for the current observation
  238 + log_header3(TERSE, get_obs_header(m_obs[i]));
271 239  
272 240 // Get CTA observation
273 241 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
274 242  
275 243 // Skip observation if it's not CTA
276 244 if (obs == NULL) {
277   - if (logTerse()) {
278   - log << " Skipping ";
279   - log << m_obs[i]->instrument();
280   - log << " observation" << std::endl;
281   - }
  245 + std::string msg = " Skipping "+m_obs[i]->instrument()+
  246 + " observation";
  247 + log_string(NORMAL, msg);
282 248 continue;
283 249 }
284 250  
... ... @@ -291,11 +257,9 @@ void ctmodel::run(void)
291 257  
292 258 // Skip observation if we have a binned observation
293 259 if (obs->eventtype() == "CountsCube") {
294   - if (logTerse()) {
295   - log << " Skipping binned ";
296   - log << obs->instrument();
297   - log << " observation" << std::endl;
298   - }
  260 + std::string msg = " Skipping binned "+m_obs[i]->instrument()+
  261 + " observation";
  262 + log_string(NORMAL, msg);
299 263 continue;
300 264 }
301 265  
... ... @@ -309,20 +273,12 @@ void ctmodel::run(void)
309 273  
310 274 } // endfor: looped over observations
311 275  
312   - // Log cube
313   - if (logTerse()) {
314   - log << std::endl;
315   - log.header1("Model cube");
316   - log << m_cube << std::endl;
317   - }
  276 + // Write model cube into header
  277 + log_header1(NORMAL, "Model cube");
  278 + log_string(NORMAL, m_cube.print());
318 279  
319   - // Restore energy dispersion flag for all CTA observations
320   - for (int i = 0; i < m_obs.size(); ++i) {
321   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
322   - if (obs != NULL) {
323   - obs->response()->apply_edisp(save_edisp[i]);
324   - }
325   - }
  280 + // Restore energy dispersion flags of all CTA observations
  281 + restore_edisp(m_obs, save_edisp);
326 282  
327 283 // Optionally publish model cube
328 284 if (m_publish) {
... ... @@ -343,27 +299,22 @@ void ctmodel::run(void)
343 299 void ctmodel::save(void)
344 300 {
345 301 // Write header
346   - if (logTerse()) {
347   - log << std::endl;
348   - log.header1("Save cube");
349   - }
  302 + log_header1(TERSE, "Save model cube");
350 303  
351 304 // Get model cube filename
352 305 m_outcube = (*this)["outcube"].filename();
353 306  
354   - // Save only if filename is non-empty
355   - if (!m_outcube.is_empty()) {
356   -
357   - // Log filename
358   - if (logTerse()) {
359   - log << gammalib::parformat("Model cube file");
360   - log << m_outcube.url() << std::endl;
361   - }
362   -
363   - // Save model cube into FITS file
  307 + // Save only if filename is non-empty and cube has some size
  308 + if (!m_outcube.is_empty() && m_cube.size() > 0) {
364 309 m_cube.save(m_outcube, clobber());
  310 + }
365 311  
  312 + // Write into logger what has been done
  313 + std::string fname = (m_outcube.is_empty()) ? "NONE" : m_outcube.url();
  314 + if (m_cube.size() == 0) {
  315 + fname.append(" (cube is empty, no file created)");
366 316 }
  317 + log_value(NORMAL, "Model cube file", fname);
367 318  
368 319 // Return
369 320 return;
... ... @@ -377,11 +328,8 @@ void ctmodel::save(void)
377 328 ***************************************************************************/
378 329 void ctmodel::publish(const std::string& name)
379 330 {
380   - // Write header
381   - if (logTerse()) {
382   - log << std::endl;
383   - log.header1("Publish model cube");
384   - }
  331 + // Write header into logger
  332 + log_header1(TERSE, "Publish model cube");
385 333  
386 334 // Set default name is user name is empty
387 335 std::string user_name(name);
... ... @@ -390,10 +338,7 @@ void ctmodel::publish(const std::string&amp; name)
390 338 }
391 339  
392 340 // Log filename
393   - if (logTerse()) {
394   - log << gammalib::parformat("Publish model cube");
395   - log << user_name << std::endl;
396   - }
  341 + log_value(NORMAL, "Model cube name", user_name);
397 342  
398 343 // Publish model cube
399 344 m_cube.counts().publish(user_name);
... ... @@ -443,6 +388,7 @@ void ctmodel::init_members(void)
443 388 m_outcube.clear();
444 389 m_apply_edisp = false;
445 390 m_publish = false;
  391 + m_chatter = static_cast<GChatter>(2);
446 392  
447 393 // Initialise protected members
448 394 m_obs.clear();
... ... @@ -468,6 +414,7 @@ void ctmodel::copy_members(const ctmodel&amp; app)
468 414 m_outcube = app.m_outcube;
469 415 m_apply_edisp = app.m_apply_edisp;
470 416 m_publish = app.m_publish;
  417 + m_chatter = app.m_chatter;
471 418  
472 419 // Copy protected members
473 420 m_obs = app.m_obs;
... ... @@ -525,7 +472,8 @@ void ctmodel::get_parameters(void)
525 472 if (obs->eventtype() == "CountsCube") {
526 473  
527 474 // Set cube from binned observation
528   - GCTAEventCube* evtcube = dynamic_cast<GCTAEventCube*>(const_cast<GEvents*>(obs->events()));
  475 + GCTAEventCube* evtcube = dynamic_cast<GCTAEventCube*>
  476 + (const_cast<GEvents*>(obs->events()));
529 477  
530 478 cube(*evtcube);
531 479  
... ... @@ -593,6 +541,7 @@ void ctmodel::get_parameters(void)
593 541  
594 542 // Get remaining parameters
595 543 m_publish = (*this)["publish"].boolean();
  544 + m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
596 545  
597 546 // Read optionally output cube filenames
598 547 if (read_ahead()) {
... ... @@ -776,97 +725,85 @@ void ctmodel::get_obs(void)
776 725 ***************************************************************************/
777 726 void ctmodel::fill_cube(const GCTAObservation* obs)
778 727 {
779   - // Continue only if observation pointer is valid
780   - if (obs != NULL) {
781   -
782   - // Get GTI and energy boundaries references for observation
783   - const GGti& gti = obs->events()->gti();
784   - const GEbounds& obs_ebounds = obs->ebounds();
785   -
786   - // Get cube energy boundaries
787   - const GEbounds& cube_ebounds = m_cube.ebounds();
  728 + // Get references to GTI and energy boundaries for the event list
  729 + const GGti& gti = obs->events()->gti();
  730 + const GEbounds& obs_ebounds = obs->ebounds();
788 731  
789   - // Initialise empty, invalid RoI
790   - GCTARoi roi;
  732 + // Get cube energy boundaries
  733 + const GEbounds& cube_ebounds = m_cube.ebounds();
791 734  
792   - // Retrieve RoI in case we have an unbinned observation
793   - if (obs->eventtype() == "EventList") {
794   - roi = obs->roi();
795   - }
  735 + // Get counts cube usage flags
  736 + std::vector<bool> usage = cube_layer_usage(cube_ebounds, obs_ebounds);
796 737  
797   - // Initialise statistics
798   - double sum = 0.0;
799   - int num_outside_ebds = 0;
800   - int num_outside_roi = 0;
  738 + // Initialise empty, invalid RoI
  739 + GCTARoi roi;
801 740  
802   - // Setup cube GTIs for this observation
803   - m_cube.gti(obs->events()->gti());
  741 + // Retrieve RoI in case we have an unbinned observation
  742 + if (obs->eventtype() == "EventList") {
  743 + roi = obs->roi();
  744 + }
804 745  
805   - // Loop over all cube bins
806   - for (int i = 0; i < m_cube.size(); ++i) {
  746 + // Initialise statistics
  747 + double sum = 0.0;
  748 + int num_outside_ebds = 0;
  749 + int num_outside_roi = 0;
807 750  
808   - // Get cube bin
809   - GCTAEventBin* bin = m_cube[i];
  751 + // Setup cube GTIs for this observation
  752 + m_cube.gti(obs->events()->gti());
810 753  
811   - // Skip bin if it is outside the energy range of the observation
812   - int index = cube_ebounds.index(bin->energy());
813   - if (index == -1 ||
814   - !obs_ebounds.contains(cube_ebounds.emin(index)+g_energy_margin,
815   - cube_ebounds.emax(index)-g_energy_margin)) {
816   - num_outside_ebds++;
817   - continue;
818   - }
  754 + // Loop over all cube bins
  755 + for (int i = 0; i < m_cube.size(); ++i) {
819 756  
820   - // Check if RoI is valid, i.e. check if we have an unbinned
821   - // observation
822   - if (roi.is_valid()) {
  757 + // Get cube bin
  758 + GCTAEventBin* bin = m_cube[i];
823 759  
824   - // Skip bin if it is outside the RoI of the observation
825   - if (!roi.contains(*bin)) {
826   - num_outside_roi++;
827   - continue;
828   - }
  760 + // Determine counts cube energy bin
  761 + int iebin = cube_ebounds.index(bin->energy());
829 762  
830   - } // endif: RoI was not valid
  763 + // Skip bin if the corresponding counts cube energy bin is not fully
  764 + // contained in the event list energy range. This avoids having
  765 + // partially filled bins.
  766 + if (!usage[iebin]) {
  767 + num_outside_ebds++;
  768 + continue;
  769 + }
831 770  
832   - // Get actual bin value
833   - double value = bin->counts();
834   -
835   - // Compute model value for cube bin
836   - double model = m_obs.models().eval(*bin, *obs) * bin->size();
  771 + // If RoI is valid then skip bin if it is outside the RoI of the
  772 + // observation
  773 + if (roi.is_valid() && !roi.contains(*bin)) {
  774 + num_outside_roi++;
  775 + continue;
  776 + }
837 777  
838   - // Add model to actual value
839   - value += model;
840   - sum += model;
  778 + // Get actual bin value
  779 + double value = bin->counts();
  780 +
  781 + // Compute model value for cube bin
  782 + double model = m_obs.models().eval(*bin, *obs) * bin->size();
841 783  
842   - // Store value
843   - bin->counts(value);
  784 + // Add model to actual value
  785 + value += model;
  786 + sum += model;
844 787  
845   - } // endfor: looped over all cube bins
  788 + // Store value
  789 + bin->counts(value);
846 790  
847   - // Append GTIs of observation to list of GTIs
848   - m_gti.extend(gti);
  791 + } // endfor: looped over all cube bins
849 792  
850   - // Update GTIs
851   - m_cube.gti(m_gti);
  793 + // Append GTIs of observation to list of GTIs
  794 + m_gti.extend(gti);
852 795  
853   - // Log results
854   - if (logTerse()) {
855   - log << gammalib::parformat("Model events in cube");
856   - log << sum << std::endl;
857   - log << gammalib::parformat("Bins outside energy range");
858   - log << num_outside_ebds << std::endl;
859   - log << gammalib::parformat("Bins outside RoI");
860   - log << num_outside_roi << std::endl;
861   - }
  796 + // Update GTIs
  797 + m_cube.gti(m_gti);
862 798  
863   - // Log cube
864   - if (logExplicit()) {
865   - log.header2("Model cube");
866   - log << m_cube << std::endl;
867   - }
  799 + // Log filling results
  800 + log_value(NORMAL, "Model events in cube", sum);
  801 + log_value(NORMAL, "Bins outside energy range", num_outside_ebds);
  802 + log_value(NORMAL, "Bins outside RoI", num_outside_roi);
868 803  
869   - } // endif: observation was valid
  804 + // Write model cube into header
  805 + log_header2(EXPLICIT, "Model cube");
  806 + log_string(EXPLICIT, m_cube.print(m_chatter));
870 807  
871 808 // Return
872 809 return;
... ...
src/ctmodel/ctmodel.hpp
  View file @0622a2c
... ... @@ -88,6 +88,7 @@ protected:
88 88 GFilename m_outcube; //!< Output model cube
89 89 bool m_apply_edisp; //!< Apply energy dispersion?
90 90 bool m_publish; //!< Publish model cube?
  91 + GChatter m_chatter; //!< Chattiness
91 92  
92 93 // Protected members
93 94 GObservations m_obs; //!< Observation container
... ...
src/ctobssim/ctobssim.cpp
  View file @0622a2c
... ... @@ -242,16 +242,15 @@ void ctobssim::run(void)
242 242 m_save_and_dispose = true;
243 243 }
244 244  
245   - // Determine the number of valid CTA observations, set energy dispersion
246   - // flag for all CTA observations and save old values in save_edisp vector
247   - int n_observations = 0;
248   - std::vector<bool> save_edisp;
249   - save_edisp.assign(m_obs.size(), false);
  245 + // Set energy dispersion flags of all CTA observations and save old
  246 + // values in save_edisp vector
  247 + std::vector<bool> save_edisp = set_edisp(m_obs, m_apply_edisp);
  248 +
  249 + // Determine the number of valid CTA observations
  250 + int n_observations = 0;
250 251 for (int i = 0; i < m_obs.size(); ++i) {
251 252 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
252 253 if (obs != NULL) {
253   - save_edisp[i] = obs->response()->apply_edisp();
254   - obs->response()->apply_edisp(m_apply_edisp);
255 254 n_observations++;
256 255 }
257 256 }
... ... @@ -262,48 +261,29 @@ void ctobssim::run(void)
262 261 m_use_xml = true;
263 262 }
264 263  
265   - // Write execution mode into logger
266   - if (logTerse()) {
267   - log << std::endl;
268   - log.header1("Execution mode");
269   - log << gammalib::parformat("Event list management");
270   - if (m_save_and_dispose) {
271   - log << "Save and dispose (reduces memory needs)" << std::endl;
272   - }
273   - else {
274   - log << "Keep events in memory" << std::endl;
275   - }
276   - log << gammalib::parformat("Output format");
277   - if (m_use_xml) {
278   - log << "Write Observation Definition XML file" << std::endl;
279   - }
280   - else {
281   - log << "Write single event list FITS file" << std::endl;
282   - }
283   - }
  264 + // Write execution model into logger
  265 + log_header1(NORMAL, "Execution mode");
  266 + std::string mode = (m_save_and_dispose)
  267 + ? "Save and dispose (reduces memory needs)"
  268 + : "Keep events in memory";
  269 + std::string xml = (m_use_xml)
  270 + ? "Write Observation Definition XML file"
  271 + : "Write single event list FITS file";
  272 + log_value(NORMAL, "Event list management", mode);
  273 + log_value(NORMAL, "Output format", xml);
284 274  
285 275 // Write seed values into logger
286   - if (logTerse()) {
287   - log << std::endl;
288   - log.header1("Seed values");
289   - for (int i = 0; i < m_rans.size(); ++i) {
290   - log << gammalib::parformat("Seed "+gammalib::str(i));
291   - log << gammalib::str(m_rans[i].seed()) << std::endl;
292   - }
  276 + log_header1(NORMAL, "Seed values");
  277 + for (int i = 0; i < m_rans.size(); ++i) {
  278 + log_value(NORMAL, "Seed "+gammalib::str(i),
  279 + gammalib::str(m_rans[i].seed()));
293 280 }
294 281  
295   - // Write observation(s) into logger
296   - if (logTerse()) {
297   - log << std::endl;
298   - log.header1(gammalib::number("Observation", m_obs.size()));
299   - log << m_obs << std::endl;
300   - }
  282 + // Write input observation container into logger
  283 + log_observations(NORMAL, m_obs, "Input observation");
301 284  
302 285 // Write header
303   - if (logTerse()) {
304   - log << std::endl;
305   - log.header1(gammalib::number("Simulate observation", m_obs.size()));
306   - }
  286 + log_header1(TERSE, gammalib::number("Simulate observation", m_obs.size()));
307 287  
308 288 // From here on the code can be parallelized if OpenMP support
309 289 // is enabled. The code in the following block corresponds to the
... ... @@ -433,13 +413,8 @@ void ctobssim::run(void)
433 413  
434 414 } // end pragma omp parallel
435 415  
436   - // Restore energy dispersion flag for all CTA observations
437   - for (int i = 0; i < m_obs.size(); ++i) {
438   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
439   - if (obs != NULL) {
440   - obs->response()->apply_edisp(save_edisp[i]);
441   - }
442   - }
  416 + // Restore energy dispersion flags of all CTA observations
  417 + restore_edisp(m_obs, save_edisp);
443 418  
444 419 // Return
445 420 return;
... ... @@ -467,16 +442,8 @@ void ctobssim::run(void)
467 442 ***************************************************************************/
468 443 void ctobssim::save(void)
469 444 {
470   - // Write header
471   - if (logTerse()) {
472   - log << std::endl;
473   - if (m_obs.size() > 1) {
474   - log.header1("Save observations");
475   - }
476   - else {
477   - log.header1("Save observation");
478   - }
479   - }
  445 + // Write header into logger
  446 + log_header1(TERSE, gammalib::number("Save observation", m_obs.size()));
480 447  
481 448 // Case A: Save event file(s) and XML metadata information
482 449 if (m_use_xml) {
... ... @@ -1460,10 +1427,7 @@ void ctobssim::save_fits(void)
1460 1427 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
1461 1428  
1462 1429 // Log filename
1463   - if (logTerse()) {
1464   - log << gammalib::parformat("Event list file");
1465   - log << m_outevents << std::endl;
1466   - }
  1430 + log_value(NORMAL, "Event list file", m_outevents);
1467 1431  
1468 1432 // Save observation into FITS file
1469 1433 obs->save(m_outevents, clobber());
... ... @@ -1492,16 +1456,7 @@ void ctobssim::save_xml(void)
1492 1456 m_prefix = (*this)["prefix"].string();
1493 1457  
1494 1458 // Issue warning if output filename has no .xml suffix
1495   - std::string suffix = gammalib::tolower(m_outevents.substr(m_outevents.length()-4,4));
1496   - if (suffix != ".xml") {
1497   - log << "*** WARNING: Name of observation definition output file \""+
1498   - m_outevents+"\"" << std::endl;
1499   - log << "*** WARNING: does not terminate with \".xml\"." << std::endl;
1500   - log << "*** WARNING: This is not an error, but might be misleading."
1501   - " It is recommended" << std::endl;
1502   - log << "*** WARNING: to use the suffix \".xml\" for observation"
1503   - " definition files." << std::endl;
1504   - }
  1459 + log_string(TERSE, warn_xml_suffix(m_outevents));
1505 1460  
1506 1461 // Save only if event lists have not yet been saved and disposed
1507 1462 if (!m_save_and_dispose) {
... ... @@ -1526,10 +1481,7 @@ void ctobssim::save_xml(void)
1526 1481 obs->eventfile(outfile);
1527 1482  
1528 1483 // Log filename
1529   - if (logTerse()) {
1530   - log << gammalib::parformat("Event list file");
1531   - log << outfile << std::endl;
1532   - }
  1484 + log_value(NORMAL, "Event list file", outfile);
1533 1485  
1534 1486 // Save observation into FITS file
1535 1487 obs->save(outfile, clobber());
... ...
src/ctselect/ctselect.cpp
  View file @0622a2c
... ... @@ -224,18 +224,11 @@ void ctselect::run(void)
224 224 // Get parameters
225 225 get_parameters();
226 226  
227   - // Write observation(s) into logger
228   - if (logTerse()) {
229   - log << std::endl;
230   - log.header1("Observations before selection");
231   - log << m_obs << std::endl;
232   - }
  227 + // Write input observation container into logger
  228 + log_observations(NORMAL, m_obs, "Input observation");
233 229  
234   - // Write header
235   - if (logTerse()) {
236   - log << std::endl;
237   - log.header1("Event selection");
238   - }
  230 + // Write header into logger
  231 + log_header1(TERSE, "Event selection");
239 232  
240 233 // Initialise counters
241 234 int n_observations = 0;
... ... @@ -243,17 +236,8 @@ void ctselect::run(void)
243 236 // Loop over all observation in the container
244 237 for (int i = 0; i < m_obs.size(); ++i) {
245 238  
246   - // Write header for observation
247   - if (logTerse()) {
248   - std::string header = m_obs[i]->instrument() + " observation";
249   - if (m_obs[i]->name().length() > 1) {
250   - header += " \"" + m_obs[i]->name() + "\"";
251   - }
252   - if (m_obs[i]->id().length() > 1) {
253   - header += " (id=" + m_obs[i]->id() +")";
254   - }
255   - log.header3(header);
256   - }
  239 + // Write header for the current observation
  240 + log_header3(TERSE, get_obs_header(m_obs[i]));
257 241  
258 242 // Initialise event input and output filenames and the event
259 243 // and GTI extension names
... ... @@ -266,21 +250,17 @@ void ctselect::run(void)
266 250  
267 251 // Skip observation if it's not CTA
268 252 if (obs == NULL) {
269   - if (logTerse()) {
270   - log << " Skipping ";
271   - log << m_obs[i]->instrument();
272   - log << " observation" << std::endl;
273   - }
  253 + std::string msg = " Skipping "+m_obs[i]->instrument()+
  254 + " observation";
  255 + log_string(NORMAL, msg);
274 256 continue;
275 257 }
276 258  
277 259 // Skip observation if we have a binned observation
278 260 if (obs->eventtype() == "CountsCube") {
279   - if (logTerse()) {
280   - log << " Skipping binned ";
281   - log << obs->instrument();
282   - log << " observation" << std::endl;
283   - }
  261 + std::string msg = " Skipping binned "+m_obs[i]->instrument()+
  262 + " observation";
  263 + log_string(NORMAL, msg);
284 264 continue;
285 265 }
286 266  
... ... @@ -297,23 +277,15 @@ void ctselect::run(void)
297 277 }
298 278 m_gtiname[i] = get_gtiname(fname.url(), m_evtname[i]);
299 279  
300   - // Log input file information
301   - if (logTerse()) {
302   - log << gammalib::parformat("Input filename");
303   - log << m_infiles[i] << std::endl;
304   - log << gammalib::parformat("Event extension name");
305   - log << m_evtname[i] << std::endl;
306   - log << gammalib::parformat("GTI extension name");
307   - log << m_gtiname[i] << std::endl;
308   - }
  280 + // Write input file information into logger
  281 + log_value(NORMAL, "Input filename", m_infiles[i]);
  282 + log_value(NORMAL, "Event extension name", m_evtname[i]);
  283 + log_value(NORMAL, "GTI extension name", m_gtiname[i]);
309 284  
310 285 // Fall through in case that the event file is empty
311 286 if (obs->events()->size() == 0) {
312   - if (logTerse()) {
313   - log << " Warning: No events in event file \"";
314   - log << m_infiles[i] << "\". Event selection skipped.";
315   - log << std::endl;
316   - }
  287 + log_string(NORMAL, " Warning: No events in event file \""+
  288 + m_infiles[i]+"\". Event selection skipped.");
317 289 continue;
318 290 }
319 291  
... ... @@ -377,11 +349,7 @@ void ctselect::run(void)
377 349 }
378 350  
379 351 // Write observation(s) into logger
380   - if (logTerse()) {
381   - log << std::endl;
382   - log.header1("Observations after selection");
383   - log << m_obs << std::endl;
384   - }
  352 + log_observations(NORMAL, m_obs, "Output observation");
385 353  
386 354 // Return
387 355 return;
... ... @@ -409,16 +377,8 @@ void ctselect::run(void)
409 377 ***************************************************************************/
410 378 void ctselect::save(void)
411 379 {
412   - // Write header
413   - if (logTerse()) {
414   - log << std::endl;
415   - if (m_obs.size() > 1) {
416   - log.header1("Save event lists");
417   - }
418   - else {
419   - log.header1("Save event list");
420   - }
421   - }
  380 + // Write header into logger
  381 + log_header1(TERSE, gammalib::number("Save event list", m_obs.size()));
422 382  
423 383 // Case A: Save event file(s) and XML metadata information
424 384 if (m_use_xml) {
... ... @@ -459,6 +419,7 @@ void ctselect::init_members(void)
459 419 m_emax = 0.0;
460 420 m_expr.clear();
461 421 m_usethres.clear();
  422 + m_chatter = static_cast<GChatter>(2);
462 423  
463 424 // Initialise protected members
464 425 m_obs.clear();
... ... @@ -496,6 +457,7 @@ void ctselect::copy_members(const ctselect&amp; app)
496 457 m_emax = app.m_emax;
497 458 m_expr = app.m_expr;
498 459 m_usethres = app.m_usethres;
  460 + m_chatter = app.m_chatter;
499 461  
500 462 // Copy protected members
501 463 m_obs = app.m_obs;
... ... @@ -610,6 +572,7 @@ void ctselect::get_parameters(void)
610 572 // Get other User parameters
611 573 m_expr = (*this)["expr"].string();
612 574 m_usethres = (*this)["usethres"].string();
  575 + m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
613 576  
614 577 // Optionally read ahead parameters so that they get correctly
615 578 // dumped into the log file
... ... @@ -726,10 +689,8 @@ void ctselect::select_events(GCTAObservation* obs,
726 689 sprintf(cmin, "%.8f", tmin);
727 690 sprintf(cmax, "%.8f", tmax);
728 691 selection = "TIME >= "+std::string(cmin)+" && TIME <= "+std::string(cmax);
729   - if (logTerse()) {
730   - log << gammalib::parformat("Time range");
731   - log << tmin << " - " << tmax << " s" << std::endl;
732   - }
  692 + log_value(NORMAL, "Time range",
  693 + gammalib::str(tmin)+" - "+gammalib::str(tmax)+" s");
733 694  
734 695 } // endif: made time selection
735 696  
... ... @@ -742,7 +703,7 @@ void ctselect::select_events(GCTAObservation* obs,
742 703 }
743 704  
744 705 // Log the requested energy selection
745   - if (logTerse()) {
  706 + if (logNormal()) {
746 707 log << gammalib::parformat("Selected energy range");
747 708 if (emin > 0.0 && emax > 0.0) {
748 709 if (emin >= emax) {
... ... @@ -804,23 +765,20 @@ void ctselect::select_events(GCTAObservation* obs,
804 765 }
805 766  
806 767 // Log the requested ROI
807   - if (logTerse()) {
808   - log << gammalib::parformat("Requested ROI");
809   - log << "Centre(RA,DEC)=(" << ra << ", " << dec << ") deg, ";
810   - log << "Radius=" << m_rad << " deg" << std::endl;
811   - }
  768 + log_value(NORMAL, "Requested RoI",
  769 + "Centre(RA,DEC)=("+gammalib::str(ra)+", "+
  770 + gammalib::str(dec)+") deg, Radius="+gammalib::str(m_rad)+
  771 + " deg");
812 772  
813 773 // If we have already an ROI then make sure that the selected
814 774 // ROI overlaps with the existing ROI
815 775 double roi_radius = list->roi().radius();
816 776 if (roi_radius > 0.0) {
817 777 GSkyDir roi_centre = list->roi().centre().dir();
818   - if (logTerse()) {
819   - log << gammalib::parformat("ROI of data");
820   - log << "Centre(RA,DEC)=(" << roi_centre.ra_deg() << ", ";
821   - log << roi_centre.dec_deg() << ") deg, ";
822   - log << "Radius=" << roi_radius << " deg" << std::endl;
823   - }
  778 + log_value(NORMAL, "RoI of data",
  779 + "Centre(RA,DEC)=("+gammalib::str(roi_centre.ra_deg())+
  780 + ", "+gammalib::str(roi_centre.dec_deg())+") deg, Radius="+
  781 + gammalib::str(roi_radius)+" deg");
824 782 GSkyDir centre;
825 783 centre.radec_deg(ra, dec);
826 784 double distance = centre.dist_deg(roi_centre);
... ... @@ -830,8 +788,8 @@ void ctselect::select_events(GCTAObservation* obs,
830 788 }
831 789  
832 790 // Log the selected ROI
833   - if (logTerse()) {
834   - log << gammalib::parformat("Selected ROI");
  791 + if (logNormal()) {
  792 + log << gammalib::parformat("Selected RoI");
835 793 if (rad <= 0.0) {
836 794 log << "None. There is no overlap between existing ";
837 795 log << "and requested ROI." << std::endl;
... ... @@ -869,10 +827,7 @@ void ctselect::select_events(GCTAObservation* obs,
869 827 } // endif: made expression selection
870 828  
871 829 // Dump cfitsio selection string
872   - if (logTerse()) {
873   - log << gammalib::parformat("cfitsio selection");
874   - log << selection << std::endl;
875   - }
  830 + log_value(NORMAL, "cfitsio selection", selection);
876 831  
877 832 // Build input filename including selection expression
878 833 std::string expression = filename;
... ... @@ -881,20 +836,14 @@ void ctselect::select_events(GCTAObservation* obs,
881 836 }
882 837  
883 838 // Dump FITS filename including selection expression
884   - if (logTerse()) {
885   - log << gammalib::parformat("FITS filename");
886   - log << expression << std::endl;
887   - }
  839 + log_value(NORMAL, "FITS filename", expression);
888 840  
889 841 // Open FITS file
890 842 GFits file(expression);
891 843  
892 844 // Log selected FITS file
893   - if (logExplicit()) {
894   - log << std::endl;
895   - log.header1("FITS file content after selection");
896   - log << file << std::endl;
897   - }
  845 + log_header1(EXPLICIT, "FITS file content after selection");
  846 + log_string(EXPLICIT, file.print(m_chatter));
898 847  
899 848 // Check if we have an events HDU
900 849 if (!file.contains(evtname)) {
... ... @@ -1022,7 +971,8 @@ GEbounds ctselect::set_ebounds(GCTAObservation* obs, const GEbounds&amp; ebounds) co
1022 971  
1023 972 // Get CTA IRF respsonse pointer and throw exception if we don't
1024 973 // find IRF
1025   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs->response());
  974 + const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
  975 + (obs->response());
1026 976 if (rsp == NULL) {
1027 977 std::string msg = "No IRF response attached to given observation.";
1028 978 throw GException::invalid_value(G_SET_EBOUNDS, msg);
... ... @@ -1268,10 +1218,7 @@ void ctselect::save_fits(void)
1268 1218 }
1269 1219  
1270 1220 // Log filename
1271   - if (logTerse()) {
1272   - log << gammalib::parformat("Event list file");
1273   - log << outfile << std::endl;
1274   - }
  1221 + log_value(NORMAL, "Event list file", outfile);
1275 1222  
1276 1223 // Save event list
1277 1224 save_event_list(obs, m_infiles[0], m_evtname[0], m_gtiname[0],
... ... @@ -1300,19 +1247,10 @@ void ctselect::save_xml(void)
1300 1247 {
1301 1248 // Get output filename and prefix
1302 1249 m_outobs = (*this)["outobs"].filename();
1303   - m_prefix = (*this)["prefix"].string();
  1250 + m_prefix = (*this)["prefix"].string();
1304 1251  
1305 1252 // Issue warning if output filename has no .xml suffix
1306   - std::string suffix = gammalib::tolower(m_outobs.substr(m_outobs.length()-4,4));
1307   - if (suffix != ".xml") {
1308   - log << "*** WARNING: Name of observation definition output file \""+
1309   - m_outobs+"\"" << std::endl;
1310   - log << "*** WARNING: does not terminate with \".xml\"." << std::endl;
1311   - log << "*** WARNING: This is not an error, but might be misleading."
1312   - " It is recommended" << std::endl;
1313   - log << "*** WARNING: to use the suffix \".xml\" for observation"
1314   - " definition files." << std::endl;
1315   - }
  1253 + log_string(TERSE, warn_xml_suffix(m_outobs));
1316 1254  
1317 1255 // Loop over all observation in the container
1318 1256 for (int i = 0; i < m_obs.size(); ++i) {
... ... @@ -1337,10 +1275,7 @@ void ctselect::save_xml(void)
1337 1275 outfile += "["+m_evtname[i]+"]";
1338 1276  
1339 1277 // Log filename
1340   - if (logTerse()) {
1341   - log << gammalib::parformat("Event list file");
1342   - log << outfile << std::endl;
1343   - }
  1278 + log_value(NORMAL, "Event list file", outfile);
1344 1279  
1345 1280 // Store output file name in observation
1346 1281 obs->eventfile(outfile);
... ... @@ -1351,10 +1286,8 @@ void ctselect::save_xml(void)
1351 1286  
1352 1287 } // endfor: looped over observations
1353 1288  
1354   - // Dump filename
1355   - if (logTerse()) {
1356   - log << "Save \""+m_outobs+"\"" << std::endl;
1357   - }
  1289 + // Write observation definition XML file name into logger
  1290 + log_value(NORMAL, "Obs. definition file", m_outobs);
1358 1291  
1359 1292 // Save observations in XML file
1360 1293 m_obs.save(m_outobs);
... ... @@ -1421,7 +1354,8 @@ void ctselect::save_event_list(const GCTAObservation* obs,
1421 1354 std::string outevt = evtname;
1422 1355 std::string outgti = gtiname;
1423 1356 if (outname.has_extname()) {
1424   - std::vector<std::string> extnames = gammalib::split(outname.extname(), ";");
  1357 + std::vector<std::string> extnames =
  1358 + gammalib::split(outname.extname(), ";");
1425 1359 if (extnames.size() > 0) {
1426 1360 std::string extname = gammalib::strip_whitespace(extnames[0]);
1427 1361 if (!extname.empty()) {
... ...
src/ctselect/ctselect.hpp
  View file @0622a2c
... ... @@ -36,7 +36,7 @@
36 36  
37 37 /* __Definitions _________________________________________________________ */
38 38 #define CTSELECT_NAME "ctselect"
39   -#define CTSELECT_VERSION "1.0.0"
  39 +#define CTSELECT_VERSION "1.2.0"
40 40  
41 41  
42 42 /***********************************************************************//**
... ... @@ -102,6 +102,7 @@ protected:
102 102 double m_emax; //!< Upper energy
103 103 std::string m_expr; //!< Selection expression
104 104 std::string m_usethres; //!< Energy threshold type
  105 + GChatter m_chatter; //!< Chattiness
105 106  
106 107 // Protected members
107 108 GObservations m_obs; //!< Observations container
... ...
src/ctskymap/ctskymap.cpp
  View file @0622a2c
... ... @@ -221,64 +221,34 @@ void ctskymap::run(void)
221 221 // Get parameters
222 222 get_parameters();
223 223  
224   - // Write observation(s) into logger
225   - if (logTerse()) {
226   - log << std::endl;
227   - if (m_obs.size() > 1) {
228   - log.header1("Observations");
229   - }
230   - else {
231   - log.header1("Observation");
232   - }
233   - log << m_obs << std::endl;
234   - }
  224 + // Write input observation container into logger
  225 + log_observations(NORMAL, m_obs, "Input observation");
235 226  
236   - // Write header
237   - if (logTerse()) {
238   - log << std::endl;
239   - if (m_obs.size() > 1) {
240   - log.header1("Map observations");
241   - }
242   - else {
243   - log.header1("Map observation");
244   - }
245   - }
  227 + // Write header into logger
  228 + log_header1(TERSE, gammalib::number("Map observation", m_obs.size()));
246 229  
247 230 // Loop over all observation in the container
248 231 for (int i = 0; i < m_obs.size(); ++i) {
249 232  
250   - // Write header for observation
251   - if (logTerse()) {
252   - std::string header = m_obs[i]->instrument() + " observation";
253   - if (m_obs[i]->name().length() > 1) {
254   - header += " \"" + m_obs[i]->name() + "\"";
255   - }
256   - if (m_obs[i]->id().length() > 1) {
257   - header += " (id=" + m_obs[i]->id() +")";
258   - }
259   - log.header3(header);
260   - }
  233 + // Write header for the current observation
  234 + log_header3(TERSE, get_obs_header(m_obs[i]));
261 235  
262 236 // Get CTA observation
263 237 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
264 238  
265 239 // Skip observation if it's not CTA
266 240 if (obs == NULL) {
267   - if (logTerse()) {
268   - log << " Skipping ";
269   - log << m_obs[i]->instrument();
270   - log << " observation" << std::endl;
271   - }
  241 + std::string msg = " Skipping "+m_obs[i]->instrument()+
  242 + " observation";
  243 + log_string(NORMAL, msg);
272 244 continue;
273 245 }
274 246  
275 247 // Skip observation if we have a binned observation
276 248 if (obs->eventtype() == "CountsCube") {
277   - if (logTerse()) {
278   - log << " Skipping binned ";
279   - log << obs->instrument();
280   - log << " observation" << std::endl;
281   - }
  249 + std::string msg = " Skipping binned "+obs->instrument()+
  250 + " observation";
  251 + log_string(NORMAL, msg);
282 252 continue;
283 253 }
284 254  
... ... @@ -309,22 +279,16 @@ void ctskymap::run(void)
309 279 void ctskymap::save(void)
310 280 {
311 281 // Write header
312   - if (logTerse()) {
313   - log << std::endl;
314   - log.header1("Save sky map");
315   - }
  282 + log_header1(TERSE, "Save sky map");
316 283  
317 284 // Get sky map filename
318 285 m_outmap = (*this)["outmap"].filename();
319 286  
320   - // Save only if filename is non-empty
321   - if (!m_outmap.is_empty()) {
  287 + // Determine whether the map is empty
  288 + bool map_is_empty = ((m_skymap.nx() == 0) || (m_skymap.ny() == 0));
322 289  
323   - // Log filename
324   - if (logTerse()) {
325   - log << gammalib::parformat("Sky map file");
326   - log << m_outmap.url() << std::endl;
327   - }
  290 + // Save sky map if filename and the map are not empty
  291 + if (!m_outmap.is_empty() && !map_is_empty) {
328 292  
329 293 // Create empty FITS file
330 294 GFits fits;
... ... @@ -337,6 +301,13 @@ void ctskymap::save(void)
337 301  
338 302 }
339 303  
  304 + // Write into logger what has been done
  305 + std::string fname = (m_outmap.is_empty()) ? "NONE" : m_outmap.url();
  306 + if (map_is_empty) {
  307 + fname.append(" (map is empty, no file created)");
  308 + }
  309 + log_value(NORMAL, "Sky map file", fname);
  310 +
340 311 // Return
341 312 return;
342 313 }
... ... @@ -349,11 +320,8 @@ void ctskymap::save(void)
349 320 ***************************************************************************/
350 321 void ctskymap::publish(const std::string& name)
351 322 {
352   - // Write header
353   - if (logTerse()) {
354   - log << std::endl;
355   - log.header1("Publish sky map");
356   - }
  323 + // Write header into logger
  324 + log_header1(TERSE, "Publish sky map");
357 325  
358 326 // Set default name is user name is empty
359 327 std::string user_name(name);
... ... @@ -361,11 +329,8 @@ void ctskymap::publish(const std::string&amp; name)
361 329 user_name = CTSKYMAP_NAME;
362 330 }
363 331  
364   - // Log filename
365   - if (logTerse()) {
366   - log << gammalib::parformat("Publish sky map");
367   - log << user_name << std::endl;
368   - }
  332 + // Write sky map name into logger
  333 + log_value(NORMAL, "Sky map name", user_name);
369 334  
370 335 // Publish sky map
371 336 m_skymap.publish(user_name);
... ... @@ -375,6 +340,60 @@ void ctskymap::publish(const std::string&amp; name)
375 340 }
376 341  
377 342  
  343 +/*==========================================================================
  344 + = =
  345 + = Private methods =
  346 + = =
  347 + ==========================================================================*/
  348 +
  349 +/***********************************************************************//**
  350 + * @brief Initialise class members
  351 + ***************************************************************************/
  352 +void ctskymap::init_members(void)
  353 +{
  354 + // Initialise members
  355 + m_obs.clear();
  356 + m_skymap.clear();
  357 + m_emin = 0.0;
  358 + m_emax = 0.0;
  359 + m_publish = false;
  360 + m_chatter = static_cast<GChatter>(2);
  361 +
  362 + // Return
  363 + return;
  364 +}
  365 +
  366 +
  367 +/***********************************************************************//**
  368 + * @brief Copy class members
  369 + *
  370 + * @param[in] app Application.
  371 + ***************************************************************************/
  372 +void ctskymap::copy_members(const ctskymap& app)
  373 +{
  374 + // Copy attributes
  375 + m_obs = app.m_obs;
  376 + m_skymap = app.m_skymap;
  377 + m_emin = app.m_emin;
  378 + m_emax = app.m_emax;
  379 + m_publish = app.m_publish;
  380 + m_chatter = app.m_chatter;
  381 +
  382 + // Return
  383 + return;
  384 +}
  385 +
  386 +
  387 +/***********************************************************************//**
  388 + * @brief Delete class members
  389 + ***************************************************************************/
  390 +void ctskymap::free_members(void)
  391 +{
  392 + // Return
  393 + return;
  394 +}
  395 +
  396 +
378 397 /***********************************************************************//**
379 398 * @brief Get application parameters
380 399 *
... ... @@ -408,6 +427,7 @@ void ctskymap::get_parameters(void)
408 427 m_emin = (*this)["emin"].real();
409 428 m_emax = (*this)["emax"].real();
410 429 m_publish = (*this)["publish"].boolean();
  430 + m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
411 431  
412 432 // Read ahead parameters
413 433 if (read_ahead()) {
... ... @@ -436,131 +456,66 @@ void ctskymap::get_parameters(void)
436 456 ***************************************************************************/
437 457 void ctskymap::map_events(GCTAObservation* obs)
438 458 {
439   - // Continue only if observation pointer is valid
440   - if (obs != NULL) {
  459 + // Make sure that the observation holds a CTA event list. If this
  460 + // is not the case then throw an exception.
  461 + if (dynamic_cast<const GCTAEventList*>(obs->events()) == NULL) {
  462 + throw GException::no_list(G_BIN_EVENTS);
  463 + }
441 464  
442   - // Make sure that the observation holds a CTA event list. If this
443   - // is not the case then throw an exception.
444   - if (dynamic_cast<const GCTAEventList*>(obs->events()) == NULL) {
445   - throw GException::no_list(G_BIN_EVENTS);
  465 + // Setup energy range covered by data
  466 + GEnergy emin;
  467 + GEnergy emax;
  468 + GEbounds ebds;
  469 + emin.TeV(m_emin);
  470 + emax.TeV(m_emax);
  471 +
  472 + // Initialise binning statistics
  473 + int num_outside_map = 0;
  474 + int num_outside_erange = 0;
  475 + int num_in_map = 0;
  476 +
  477 + // Fill sky map
  478 + GCTAEventList* events = static_cast<GCTAEventList*>
  479 + (const_cast<GEvents*>(obs->events()));
  480 + for (int i = 0; i < events->size(); ++i) {
  481 +
  482 + // Get event
  483 + GCTAEventAtom* event = (*events)[i];
  484 +
  485 + // Skip if energy is out of range
  486 + if (event->energy() < emin || event->energy() > emax) {
  487 + num_outside_erange++;
  488 + continue;
446 489 }
447 490  
448   - // Setup energy range covered by data
449   - GEnergy emin;
450   - GEnergy emax;
451   - GEbounds ebds;
452   - emin.TeV(m_emin);
453   - emax.TeV(m_emax);
454   -
455   - // Initialise binning statistics
456   - int num_outside_map = 0;
457   - int num_outside_erange = 0;
458   - int num_in_map = 0;
459   -
460   - // Fill sky map
461   - GCTAEventList* events = static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
462   - for (int i = 0; i < events->size(); ++i) {
463   -
464   - // Get event
465   - GCTAEventAtom* event = (*events)[i];
466   -
467   - // Skip if energy is out of range
468   - if (event->energy() < emin || event->energy() > emax) {
469   - num_outside_erange++;
470   - continue;
471   - }
472   -
473   - // Determine sky pixel
474   - GCTAInstDir* inst = (GCTAInstDir*)&(event->dir());
475   - GSkyDir dir = inst->dir();
476   - GSkyPixel pixel = m_skymap.dir2pix(dir);
477   -
478   - // Skip if pixel is out of range
479   - if (pixel.x() < -0.5 || pixel.x() > (m_skymap.nx() - 0.5) ||
480   - pixel.y() < -0.5 || pixel.y() > (m_skymap.ny() - 0.5)) {
481   - num_outside_map++;
482   - continue;
483   - }
484   -
485   - // Fill event in skymap
486   - m_skymap(pixel, 0) += 1.0;
487   - num_in_map++;
488   -
489   - } // endfor: looped over all events
490   -
491   - // Log binning results
492   - if (logTerse()) {
493   - log << gammalib::parformat("Events in list");
494   - log << obs->events()->size() << std::endl;
495   - log << gammalib::parformat("Events in map");
496   - log << num_in_map << std::endl;
497   - log << gammalib::parformat("Events outside map area");
498   - log << num_outside_map << std::endl;
499   - log << gammalib::parformat("Events outside energies");
500   - log << num_outside_erange << std::endl;
501   - }
  491 + // Determine sky pixel
  492 + GCTAInstDir* inst = (GCTAInstDir*)&(event->dir());
  493 + GSkyDir dir = inst->dir();
  494 + GSkyPixel pixel = m_skymap.dir2pix(dir);
502 495  
503   - // Log map
504   - if (logVerbose()) {
505   - log.indent(1);
506   - log << m_skymap << std::endl;
507   - log.indent(0);
  496 + // Skip if pixel is out of range
  497 + if (pixel.x() < -0.5 || pixel.x() > (m_skymap.nx() - 0.5) ||
  498 + pixel.y() < -0.5 || pixel.y() > (m_skymap.ny() - 0.5)) {
  499 + num_outside_map++;
  500 + continue;
508 501 }
509 502  
510   - } // endif: observation was valid
  503 + // Fill event in skymap
  504 + m_skymap(pixel, 0) += 1.0;
  505 + num_in_map++;
511 506  
512   - // Return
513   - return;
514   -}
  507 + } // endfor: looped over all events
515 508  
  509 + // Log binning results
  510 + log_value(NORMAL, "Events in list", obs->events()->size());
  511 + log_value(NORMAL, "Events in map", num_in_map);
  512 + log_value(NORMAL, "Events outside map area", num_outside_map);
  513 + log_value(NORMAL, "Events outside energies", num_outside_erange);
516 514  
517   -/*==========================================================================
518   - = =
519   - = Private methods =
520   - = =
521   - ==========================================================================*/
  515 + // Write sky map into header
  516 + log_header2(EXPLICIT, "Sky map");
  517 + log_string(EXPLICIT, m_skymap.print(m_chatter));
522 518  
523   -/***********************************************************************//**
524   - * @brief Initialise class members
525   - ***************************************************************************/
526   -void ctskymap::init_members(void)
527   -{
528   - // Initialise members
529   - m_obs.clear();
530   - m_skymap.clear();
531   - m_emin = 0.0;
532   - m_emax = 0.0;
533   - m_publish = false;
534   -
535   - // Return
536   - return;
537   -}
538   -
539   -
540   -/***********************************************************************//**
541   - * @brief Copy class members
542   - *
543   - * @param[in] app Application.
544   - ***************************************************************************/
545   -void ctskymap::copy_members(const ctskymap& app)
546   -{
547   - // Copy attributes
548   - m_obs = app.m_obs;
549   - m_skymap = app.m_skymap;
550   - m_emin = app.m_emin;
551   - m_emax = app.m_emax;
552   - m_publish = app.m_publish;
553   -
554   - // Return
555   - return;
556   -}
557   -
558   -
559   -/***********************************************************************//**
560   - * @brief Delete class members
561   - ***************************************************************************/
562   -void ctskymap::free_members(void)
563   -{
564 519 // Return
565 520 return;
566 521 }
... ...
src/ctskymap/ctskymap.hpp
  View file @0622a2c
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTSKYMAP_NAME "ctskymap"
37   -#define CTSKYMAP_VERSION "1.1.0"
  37 +#define CTSKYMAP_VERSION "1.2.0"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -77,6 +77,7 @@ protected:
77 77 double m_emin; //!< Minimum energy (TeV)
78 78 double m_emax; //!< Maximum energy (TeV)
79 79 bool m_publish; //!< Publish sky map?
  80 + GChatter m_chatter; //!< Chattiness
80 81  
81 82 // Protected members
82 83 GObservations m_obs; //!< Observation container
... ...
src/cttsmap/cttsmap.cpp
  View file @0622a2c
... ... @@ -215,36 +215,15 @@ void cttsmap::run(void)
215 215 // Get task parameters
216 216 get_parameters();
217 217  
218   - // Set energy dispersion flag for all CTA observations and save old
  218 + // Set energy dispersion flags of all CTA observations and save old
219 219 // values in save_edisp vector
220   - std::vector<bool> save_edisp;
221   - save_edisp.assign(m_obs.size(), false);
222   - for (int i = 0; i < m_obs.size(); ++i) {
223   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
224   - if (obs != NULL) {
225   - save_edisp[i] = obs->response()->apply_edisp();
226   - obs->response()->apply_edisp(m_apply_edisp);
227   - }
228   - }
  220 + std::vector<bool> save_edisp = set_edisp(m_obs, m_apply_edisp);
229 221  
230   - // Write observation(s) into logger
231   - if (logTerse()) {
232   - log << std::endl;
233   - if (m_obs.size() > 1) {
234   - log.header1("Observations");
235   - }
236   - else {
237   - log.header1("Observation");
238   - }
239   - log << m_obs << std::endl;
240   - }
  222 + // Write input observation container into logger
  223 + log_observations(NORMAL, m_obs, "Input observation");
241 224  
242 225 // Write header
243   - if (logTerse()) {
244   - log << std::endl;
245   - log.header1("Initialise TS map");
246   - }
247   -
  226 + log_header1(TERSE, "Initialise TS map");
248 227  
249 228 // Get bins to be computed
250 229 // To split the computation on several jobs,
... ... @@ -270,10 +249,7 @@ void cttsmap::run(void)
270 249 if (m_logL0 == 0.0) {
271 250  
272 251 // Write header
273   - if (logTerse()) {
274   - log << std::endl;
275   - log.header1("Compute NULL Hypothesis for TS computation");
276   - }
  252 + log_header1(TERSE, "Compute NULL Hypothesis for TS computation");
277 253  
278 254 // Compute likelihood without the test source
279 255 m_obs.models(models);
... ... @@ -283,10 +259,7 @@ void cttsmap::run(void)
283 259 }
284 260  
285 261 // Write header
286   - if (logTerse()) {
287   - log << std::endl;
288   - log.header1("Generate TS map");
289   - }
  262 + log_header1(TERSE, "Generate TS map");
290 263  
291 264 // Loop over grid positions
292 265 for (int i = binmin; i < binmax; ++i) {
... ... @@ -295,12 +268,8 @@ void cttsmap::run(void)
295 268 GSkyDir bincentre = m_tsmap.inx2dir(i);
296 269  
297 270 // Header for verbose logging
298   - if (logExplicit()) {
299   - std::string msg = "Computing TS for bin number "+gammalib::str(i)+
300   - " at "+bincentre.print();
301   - log << std::endl;
302   - log.header2(msg);
303   - }
  271 + log_header2(EXPLICIT, "Computing TS for bin number "+gammalib::str(i)+
  272 + " at "+bincentre.print());
304 273  
305 274 // Add test source at current bin position
306 275 if (m_testsource != NULL) {
... ... @@ -360,13 +329,8 @@ void cttsmap::run(void)
360 329 // Bring models to initial state
361 330 m_obs.models(models_orig);
362 331  
363   - // Restore energy dispersion flag for all CTA observations
364   - for (int i = 0; i < m_obs.size(); ++i) {
365   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
366   - if (obs != NULL) {
367   - obs->response()->apply_edisp(save_edisp[i]);
368   - }
369   - }
  332 + // Restore energy dispersion flags of all CTA observations
  333 + restore_edisp(m_obs, save_edisp);
370 334  
371 335 // Optionally publish TS map
372 336 if (m_publish) {
... ... @@ -386,22 +350,16 @@ void cttsmap::run(void)
386 350 void cttsmap::save(void)
387 351 {
388 352 // Write header
389   - if (logTerse()) {
390   - log << std::endl;
391   - log.header1("Save TS map");
392   - }
  353 + log_header1(TERSE, "Save TS map");
393 354  
394 355 // Get TS map filename
395 356 m_outmap = (*this)["outmap"].filename();
396 357  
397   - // Save only if filename is non-empty
398   - if (!m_outmap.is_empty()) {
  358 + // Determine whether the TS map is empty
  359 + bool map_is_empty = ((m_tsmap.nx() == 0) || (m_tsmap.ny() == 0));
399 360  
400   - // Log filename
401   - if (logTerse()) {
402   - log << gammalib::parformat("TS map file");
403   - log << m_outmap.url() << std::endl;
404   - }
  361 + // Save only if filename is non-empty
  362 + if (!m_outmap.is_empty() && !map_is_empty) {
405 363  
406 364 // Create fits file
407 365 GFits fits;
... ... @@ -428,6 +386,13 @@ void cttsmap::save(void)
428 386  
429 387 } // endif: TS map filename was valid
430 388  
  389 + // Write into logger what has been done
  390 + std::string fname = (m_outmap.is_empty()) ? "NONE" : m_outmap.url();
  391 + if (map_is_empty) {
  392 + fname.append(" (TS map is empty, no file created)");
  393 + }
  394 + log_value(NORMAL, "TS map file", fname);
  395 +
431 396 // Return
432 397 return;
433 398 }
... ... @@ -440,11 +405,8 @@ void cttsmap::save(void)
440 405 ***************************************************************************/
441 406 void cttsmap::publish(const std::string& name)
442 407 {
443   - // Write header
444   - if (logTerse()) {
445   - log << std::endl;
446   - log.header1("Publish TS map");
447   - }
  408 + // Write header into logger
  409 + log_header1(TERSE, "Publish TS map");
448 410  
449 411 // Set default name is user name is empty
450 412 std::string user_name(name);
... ... @@ -452,11 +414,8 @@ void cttsmap::publish(const std::string&amp; name)
452 414 user_name = CTTSMAP_NAME;
453 415 }
454 416  
455   - // Log filename
456   - if (logTerse()) {
457   - log << gammalib::parformat("Publish TS map");
458   - log << user_name << std::endl;
459   - }
  417 + // Write TS map name into logger
  418 + log_value(NORMAL, "TS map name", user_name);
460 419  
461 420 // Publish TS map
462 421 m_tsmap.publish(user_name);
... ...
src/cttsmap/cttsmap.hpp
  View file @0622a2c
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTTSMAP_NAME "cttsmap"
37   -#define CTTSMAP_VERSION "1.1.0"
  37 +#define CTTSMAP_VERSION "1.2.0"
38 38  
39 39  
40 40 /***********************************************************************//**
... ...
src/support/ctool.cpp
  View file @0622a2c
... ... @@ -1543,6 +1543,49 @@ GEnergies ctool::insert_energy_boundaries(const GEnergies&amp; energies,
1543 1543  
1544 1544  
1545 1545 /***********************************************************************//**
  1546 + * @brief Determine the counts cube layer usage.
  1547 + *
  1548 + * @param[in] ebounds_cube Energy boundaries of the counts cube.
  1549 + * @param[in] ebounds_list Energy boundaries of the event list.
  1550 + * @return Vector of usage flags.
  1551 + *
  1552 + * Determines a vector of counts cube layer usage flags that signal whether
  1553 + * an event should actually be filled in the counts cube or not. This makes
  1554 + * sure that no partially filled bins will exist in the counts cube.
  1555 + ***************************************************************************/
  1556 +std::vector<bool> ctool::cube_layer_usage(const GEbounds& cube_ebounds,
  1557 + const GEbounds& list_ebounds) const
  1558 +{
  1559 + // Set energy margin
  1560 + const GEnergy energy_margin(0.01, "GeV");
  1561 +
  1562 + // Initialise usage vector
  1563 + std::vector<bool> usage(cube_ebounds.size(), true);
  1564 +
  1565 + // Loop over all energy bins of the cube
  1566 + for (int iebin = 0; iebin < cube_ebounds.size(); ++iebin) {
  1567 +
  1568 + // If the counts cube energy bin is fully contained within the
  1569 + // energy boundaries of the event list the signal usage of this
  1570 + // counts cube bin. Partially overlapping energy bins are signaled
  1571 + // for non-usage, which avoids having partially filled bins in the
  1572 + // counts cube. Some margin is applied that effectively reduces the
  1573 + // width of the counts cube energy bin. This should cope with any
  1574 + // rounding errors that come from reading and writing the energy
  1575 + // boundary information to a FITS file.
  1576 + usage[iebin] = list_ebounds.contains(cube_ebounds.emin(iebin) +
  1577 + energy_margin,
  1578 + cube_ebounds.emax(iebin) -
  1579 + energy_margin);
  1580 +
  1581 + } // endfor: looped over all energy bins
  1582 +
  1583 + // Return usage
  1584 + return usage;
  1585 +}
  1586 +
  1587 +
  1588 +/***********************************************************************//**
1546 1589 * @brief Set warning string if there are too few energies
1547 1590 *
1548 1591 * @param[in] energies Energies.
... ... @@ -1584,6 +1627,37 @@ std::string ctool::warn_too_few_energies(const GEnergies&amp; energies) const
1584 1627 }
1585 1628  
1586 1629  
  1630 +/***********************************************************************//**
  1631 + * @brief Set warning string if file has no .xml suffix
  1632 + *
  1633 + * @param[in] filename Filename.
  1634 + * @return Warning string.
  1635 + *
  1636 + * Sets a warning string if the filename has not ".xml" suffix.
  1637 + ***************************************************************************/
  1638 +std::string ctool::warn_xml_suffix(const GFilename& filename) const
  1639 +{
  1640 + // Initialise warning
  1641 + std::string warning;
  1642 +
  1643 + // Get filename suffix
  1644 + std::string fname = std::string(filename);
  1645 + std::string suffix = gammalib::tolower(fname.substr(fname.length()-4,4));
  1646 +
  1647 + // Warn if there are not enough energy bins
  1648 + if (suffix != ".xml") {
  1649 + warning.append("\nWARNING: Name of observation definition output file "
  1650 + "\""+filename+"\" does not terminate with \".xml\"."
  1651 + "\n This is not an error, but may be misleading. "
  1652 + "It is recommended to use the suffix"
  1653 + "\n \".xml\" for observation definition files.");
  1654 + }
  1655 +
  1656 + // Return warning
  1657 + return warning;
  1658 +}
  1659 +
  1660 +
1587 1661 /*==========================================================================
1588 1662 = =
1589 1663 = Private methods for SWIG =
... ...
src/support/ctool.hpp
  View file @0622a2c
... ... @@ -105,9 +105,12 @@ protected:
105 105 std::string get_obs_header(const GObservation* obs);
106 106 GEnergies insert_energy_boundaries(const GEnergies& energies,
107 107 const GCTAObservation& obs);
  108 + std::vector<bool> cube_layer_usage(const GEbounds& cube_ebounds,
  109 + const GEbounds& list_ebounds) const;
108 110  
109 111 // Protected warning strings
110 112 std::string warn_too_few_energies(const GEnergies& energies) const;
  113 + std::string warn_xml_suffix(const GFilename& filename) const;
111 114  
112 115 // Protected members
113 116 bool m_read_ahead; //!< Read ahead output parameters
... ...
test/test_ctmapcube.py
  View file @0622a2c
... ... @@ -18,6 +18,7 @@
18 18 # along with this program. If not, see <http://www.gnu.org/licenses/>.
19 19 #
20 20 # ==========================================================================
  21 +import os
21 22 import gammalib
22 23 import ctools
23 24 from testing import test
... ... @@ -85,7 +86,7 @@ class Test(test):
85 86 'Check successful execution from command line')
86 87  
87 88 # Check map cube
88   - self._check_cube('ctmapcube_cmd1.fits')
  89 + self._check_result_file('ctmapcube_cmd1.fits')
89 90  
90 91 # Setup ctmapcube command
91 92 cmd = ctmapcube+' inmodel="events_that_do_not_exist.fits"'+ \
... ... @@ -107,10 +108,26 @@ class Test(test):
107 108 """
108 109 Test ctmapcube from Python
109 110 """
110   - # Set-up ctmapcube
  111 + # Allocate ctmapcube
111 112 mapcube = ctools.ctmapcube()
  113 +
  114 + # Check that empty ctmapcube tool holds an map cube that has no
  115 + # energy bins
  116 + self._check_cube(mapcube.mapcube(), nmaps=0, npixels=0)
  117 +
  118 + # Check that saving does not nothing
  119 + mapcube['outcube'] = 'ctmapcube_py0.fits'
  120 + mapcube['logfile'] = 'ctmapcube_py0.log'
  121 + mapcube.logFileOpen()
  122 + mapcube.save()
  123 + self.test_assert(not os.path.isfile('ctmapcube_py0.fits'),
  124 + 'Check that no map cube has been created')
  125 +
  126 + # Check that clearing does not lead to an exception or segfault
  127 + mapcube.clear()
  128 +
  129 + # Now set ctmapcube parameters
112 130 mapcube['inmodel'] = self._model
113   - mapcube['outcube'] = 'ctmapcube_py1.fits'
114 131 mapcube['ebinalg'] = 'LOG'
115 132 mapcube['emin'] = 0.1
116 133 mapcube['emax'] = 100.0
... ... @@ -122,6 +139,7 @@ class Test(test):
122 139 mapcube['proj'] = 'CAR'
123 140 mapcube['xref'] = 83.63
124 141 mapcube['yref'] = 22.01
  142 + mapcube['outcube'] = 'ctmapcube_py1.fits'
125 143 mapcube['logfile'] = 'ctmapcube_py1.log'
126 144 mapcube['chatter'] = 2
127 145  
... ... @@ -133,31 +151,76 @@ class Test(test):
133 151 mapcube.save()
134 152  
135 153 # Check map cube
136   - self._check_cube('ctmapcube_py1.fits')
  154 + self._check_result_file('ctmapcube_py1.fits')
  155 +
  156 + # Copy ctmapcube tool
  157 + cpy_mapcube = mapcube.copy()
  158 +
  159 + # Check map cube of ctmapcube copy
  160 + self._check_cube(cpy_mapcube.mapcube())
  161 +
  162 + # Execute copy of ctmapcube tool again, now with a higher chatter
  163 + # level than before
  164 + cpy_mapcube['outcube'] = 'ctmapcube_py2.fits'
  165 + cpy_mapcube['logfile'] = 'ctmapcube_py2.log'
  166 + cpy_mapcube['publish'] = True
  167 + cpy_mapcube['chatter'] = 3
  168 + cpy_mapcube.logFileOpen() # Needed to get a new log file
  169 + cpy_mapcube.execute()
  170 +
  171 + # Check result file
  172 + self._check_result_file('ctmapcube_py2.fits')
  173 +
  174 + # Now clear copy of ctmapcube tool
  175 + cpy_mapcube.clear()
  176 +
  177 + # Check that the cleared copy has also cleared the map cube
  178 + self._check_cube(cpy_mapcube.mapcube(), nmaps=0, npixels=0)
137 179  
138 180 # Return
139 181 return
140 182  
141   - # Check map cube
142   - def _check_cube(self, filename):
  183 + # Check result file
  184 + def _check_result_file(self, filename):
143 185 """
144 186 Check content of map cube
145 187  
146 188 Parameters
147 189 ----------
148 190 filename : str
149   - Map cube file name.
  191 + Map cube file name
150 192 """
151 193 # Load map cube
152 194 cube = gammalib.GModelSpatialDiffuseCube(filename)
153 195  
154   - # Test map cube
155   - #self.test_value(cube.maps(), 10, "10 maps")
156   - #self.test_value(cube.pixels(), 10000, "10000 map pixels")
157   - # The map cube is not loaded by default !!!! We should add the method so that
158   - # the attributes are known !!!!
159   - self.test_value(cube.maps(), 0, '10 maps')
160   - self.test_value(cube.pixels(), 0, '10000 map pixels')
  196 + # Check map cube
  197 + self._check_cube(cube)
  198 +
  199 + # Return
  200 + return
  201 +
  202 + # Check map cube
  203 + def _check_cube(self, cube, nmaps=10, npixels=10000):
  204 + """
  205 + Check map cube
  206 +
  207 + Parameters
  208 + ----------
  209 + cube : `~gammalib.GModelSpatialDiffuseCube`
  210 + Map cube
  211 + nmaps : int, optional
  212 + Number of maps
  213 + npixels : int, optional
  214 + Number of pixels
  215 + """
  216 + # Get energies (this forces loading in case the map cube is not
  217 + # loaded)
  218 + energies = cube.energies()
  219 + self.test_value(len(cube.energies()), nmaps, 'Check number of energies')
  220 +
  221 + # Check dimensions
  222 + self.test_value(cube.maps(), nmaps, 'Check number of maps')
  223 + self.test_value(cube.pixels(), npixels, 'Check number of pixels')
161 224  
162 225 # Return
163 226 return
... ...
test/test_ctmodel.py
  View file @0622a2c
... ... @@ -18,6 +18,7 @@
18 18 # along with this program. If not, see <http://www.gnu.org/licenses/>.
19 19 #
20 20 # ==========================================================================
  21 +import os
21 22 import gammalib
22 23 import ctools
23 24 from testing import test
... ... @@ -71,8 +72,8 @@ class Test(test):
71 72 ' inmodel="'+self._model+'"'+ \
72 73 ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
73 74 ' rad=5.0 ra=83.63 dec=22.01 tmin=0.0 tmax=1800'+ \
74   - ' emin=0.1 emax=100.0 enumbins=20 nxpix=200 nypix=200'+ \
75   - ' binsz=0.02 coordsys="CEL" proj="CAR"'+ \
  75 + ' emin=0.1 emax=100.0 enumbins=10 nxpix=40 nypix=40'+ \
  76 + ' binsz=0.1 coordsys="CEL" proj="CAR"'+ \
76 77 ' xref=83.63 yref=22.01'+ \
77 78 ' logfile="ctmodel_cmd1.log" chatter=1'
78 79  
... ... @@ -94,8 +95,8 @@ class Test(test):
94 95 ' inmodel="model_that_does_not_exist.xml"'+ \
95 96 ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
96 97 ' rad=5.0 ra=83.63 dec=22.01 tmin=0.0 tmax=1800'+ \
97   - ' emin=0.1 emax=100.0 enumbins=20 nxpix=200 nypix=200'+ \
98   - ' binsz=0.02 coordsys="CEL" proj="CAR"'+ \
  98 + ' emin=0.1 emax=100.0 enumbins=10 nxpix=40 nypix=40'+ \
  99 + ' binsz=0.1 coordsys="CEL" proj="CAR"'+ \
99 100 ' xref=83.63 yref=22.01'+ \
100 101 ' logfile="ctmodel_cmd2.log" chatter=1'
101 102  
... ... @@ -111,10 +112,25 @@ class Test(test):
111 112 """
112 113 Test ctmodel from Python
113 114 """
114   - # Set-up ctmodel from scratch
  115 + # Allocate ctmodel
115 116 model = ctools.ctmodel()
  117 +
  118 + # Check that empty ctmodel tool holds an empty model cube
  119 + self._check_cube(model.cube(), nx=0, ny=0, nebins=0)
  120 +
  121 + # Check that saving does not nothing
  122 + model['outcube'] = 'ctmodel_py0.fits'
  123 + model['logfile'] = 'ctmodel_py0.log'
  124 + model.logFileOpen()
  125 + model.save()
  126 + self.test_assert(not os.path.isfile('ctmodel_py0.fits'),
  127 + 'Check that no model cube has been created')
  128 +
  129 + # Check that clearing does not lead to an exception or segfault
  130 + model.clear()
  131 +
  132 + # Now set ctmodel parameters
116 133 model['incube'] = 'NONE'
117   - model['outcube'] = 'ctmodel_py1.fits'
118 134 model['inmodel'] = self._model
119 135 model['inobs'] = 'NONE'
120 136 model['expcube'] = 'NONE'
... ... @@ -128,15 +144,16 @@ class Test(test):
128 144 model['tmin'] = 0
129 145 model['tmax'] = 1800
130 146 model['emin'] = 0.1
131   - model['emax'] = 100
132   - model['enumbins'] = 20
133   - model['nxpix'] = 200
134   - model['nypix'] = 200
135   - model['binsz'] = 0.02
  147 + model['emax'] = 100.0
  148 + model['enumbins'] = 10
  149 + model['nxpix'] = 40
  150 + model['nypix'] = 40
  151 + model['binsz'] = 0.1
136 152 model['coordsys'] = 'CEL'
137 153 model['proj'] = 'CAR'
138 154 model['xref'] = 83.63
139 155 model['yref'] = 22.01
  156 + model['outcube'] = 'ctmodel_py1.fits'
140 157 model['logfile'] = 'ctmodel_py1.log'
141 158 model['chatter'] = 2
142 159  
... ... @@ -148,6 +165,30 @@ class Test(test):
148 165 # Check result file
149 166 self._check_result_file('ctmodel_py1.fits')
150 167  
  168 + # Copy ctmodel tool
  169 + cpy_model = model.copy()
  170 +
  171 + # Check model cube of ctmodel copy
  172 + self._check_cube(cpy_model.cube())
  173 +
  174 + # Execute copy of ctmodel tool again, now with a higher chatter
  175 + # level than before
  176 + cpy_model['outcube'] = 'ctmodel_py2.fits'
  177 + cpy_model['logfile'] = 'ctmodel_py2.log'
  178 + cpy_model['publish'] = True
  179 + cpy_model['chatter'] = 3
  180 + cpy_model.logFileOpen() # Needed to get a new log file
  181 + cpy_model.execute()
  182 +
  183 + # Check result file
  184 + self._check_result_file('ctmodel_py2.fits')
  185 +
  186 + # Now clear copy of ctmodel tool
  187 + cpy_model.clear()
  188 +
  189 + # Check that the cleared copy has also cleared the model cube
  190 + self._check_cube(cpy_model.cube(), nx=0, ny=0, nebins=0)
  191 +
151 192 # Return
152 193 return
153 194  
... ... @@ -155,20 +196,41 @@ class Test(test):
155 196 def _check_result_file(self, filename):
156 197 """
157 198 Check result file
  199 +
  200 + Parameters
  201 + ----------
  202 + filename : str
  203 + Model cube file name
158 204 """
159 205 # Open result file
160   - fits = gammalib.GFits(filename)
161   -
162   - # Get HDUs
163   - cube = fits['Primary']
164   - ebounds = fits['EBOUNDS']
165   - gti = fits['GTI']
166   -
167   - # Check dimensions
168   - self.test_value(cube.naxis(), 3, 'Check for 3 cube dimensions')
169   - self.test_value(cube.naxes(0), 200, 'Check for 200 pixels in X')
170   - self.test_value(cube.naxes(1), 200, 'Check for 200 pixels in Y')
171   - self.test_value(cube.naxes(2), 20, 'Check for 20 pixels in Z')
  206 + cube = gammalib.GCTAEventCube(filename)
  207 +
  208 + # Check cube
  209 + self._check_cube(cube)
  210 +
  211 + # Return
  212 + return
  213 +
  214 + # Check model cube
  215 + def _check_cube(self, cube, nx=40, ny=40, nebins=10):
  216 + """
  217 + Check model cube
  218 +
  219 + Parameters
  220 + ----------
  221 + cube : `~gammalib.GCTAEventCube`
  222 + Model cube
  223 + nx : int, optional
  224 + Number of pixels in X direction
  225 + ny : int, optional
  226 + Number of pixels in Y direction
  227 + nebins : int, optional
  228 + Number of energy bins
  229 + """
  230 + # Check number of bins per axis
  231 + self.test_value(cube.nx(), nx, 'Check number of X pixels')
  232 + self.test_value(cube.ny(), ny, 'Check number of Y pixels')
  233 + self.test_value(cube.ebins(), nebins, 'Check number of energy bins')
172 234  
173 235 # Return
174 236 return
... ...