Johan Bregeon / ctools

Commit 115eb8a9ed6d424e6043afcf214098b3e5cd3e57

Authored by Jürgen Knödlseder
1 parent b8e1637c

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

Inline Side-by-side
Showing 7 changed files with 459 additions and 271 deletions
src/cterror/cterror.cpp
  View file @115eb8a
... ... @@ -228,35 +228,15 @@ void cterror::run(void)
228 228 // Get task parameters
229 229 get_parameters();
230 230  
231   - // Set energy dispersion flag for all CTA observations and save old
  231 + // Set energy dispersion flags of all CTA observations and save old
232 232 // values in save_edisp vector
233   - std::vector<bool> save_edisp;
234   - save_edisp.assign(m_obs.size(), false);
235   - for (int i = 0; i < m_obs.size(); ++i) {
236   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
237   - if (obs != NULL) {
238   - save_edisp[i] = obs->response()->apply_edisp();
239   - obs->response()->apply_edisp(m_apply_edisp);
240   - }
241   - }
  233 + std::vector<bool> save_edisp = set_edisp(m_obs, m_apply_edisp);
242 234  
243   - // Write observation(s) into logger
244   - if (logTerse()) {
245   - log << std::endl;
246   - if (m_obs.size() > 1) {
247   - log.header1("Observations");
248   - }
249   - else {
250   - log.header1("Observation");
251   - }
252   - log << m_obs << std::endl;
253   - }
  235 + // Write input observation container into logger
  236 + log_observations(NORMAL, m_obs, "Input observation");
254 237  
255   - // Write header
256   - if (logTerse()) {
257   - log << std::endl;
258   - log.header1("Compute best-fit likelihood");
259   - }
  238 + // Write header into logger
  239 + log_header1(TERSE, "Compute best-fit likelihood");
260 240  
261 241 // Optimize and save best log-likelihood
262 242 m_obs.optimize(m_opt);
... ... @@ -266,13 +246,10 @@ void cterror::run(void)
266 246 // Store optimizer for later recovery
267 247 GOptimizerLM best_opt = m_opt;
268 248  
269   - // Write optimised model into logger
270   - if (logTerse()) {
271   - log << m_opt << std::endl;
272   - log << gammalib::parformat("Maximum log likelihood");
273   - log << gammalib::str(m_best_logL,3) << std::endl;
274   - log << m_obs.models() << std::endl;
275   - }
  249 + // Write optimisation results and models into logger
  250 + log_string(NORMAL, m_opt.print(m_chatter));
  251 + log_value(NORMAL, "Maximum log likelihood", gammalib::str(m_best_logL,3));
  252 + log_string(NORMAL, m_obs.models().print(m_chatter));
276 253  
277 254 // Continue only if source model exists
278 255 if (m_obs.models().contains(m_srcname)) {
... ... @@ -310,63 +287,39 @@ void cterror::run(void)
310 287 double parmax = std::min(m_model_par->factor_max(),
311 288 m_value + 10.0*m_model_par->factor_error());
312 289  
313   - // Write header
314   - if (logTerse()) {
315   - log << std::endl;
316   - log.header1("Compute error for source \""+m_srcname+"\""
317   - " parameter \""+m_model_par->name()+"\"");
318   - log << gammalib::parformat("Confidence level");
319   - log << m_confidence*100.0 << "%" << std::endl;
320   - log << gammalib::parformat("Log-likelihood difference");
321   - log << m_dlogL << std::endl;
322   - log << gammalib::parformat("Initial factor range");
323   - log << "[";
324   - log << parmin;
325   - log << ", ";
326   - log << parmax;
327   - log << "]" << std::endl;
328   - }
  290 + // Write header and initial parameters into logger
  291 + log_header1(TERSE, "Compute error for source \""+m_srcname+"\""
  292 + " parameter \""+m_model_par->name()+"\"");
  293 + log_value(NORMAL, "Confidence level",
  294 + gammalib::str(m_confidence*100.0)+" %");
  295 + log_value(NORMAL, "Log-likelihood difference", m_dlogL);
  296 + log_value(NORMAL, "Initial factor range",
  297 + "["+gammalib::str(parmin)+", "+gammalib::str(parmax)+"]");
329 298  
330   - // Compute lower boundary
  299 + // Compute lower and upper boundaries
331 300 double value_lo = error_bisection(parmin, m_value);
332   -
333   - // Write lower parameter value
334   - if (logTerse()) {
335   - log << gammalib::parformat("Lower parameter factor");
336   - log << value_lo << std::endl;
337   - }
338   -
339   - // Compute upper boundary
340 301 double value_hi = error_bisection(m_value, parmax);
341 302  
342   - // Write upper parameter value
343   - if (logTerse()) {
344   - log << gammalib::parformat("Upper parameter factor");
345   - log << value_hi << std::endl;
346   - }
347   -
348 303 // Compute errors
349   - double error = 0.5 * (value_hi - value_lo);
350   - double error_neg = m_value - value_lo;
351   - double error_pos = value_hi - m_value;
352   - //double error_max = std::max(value_hi-m_value, m_value-value_lo);
353   - //double error_min = std::min(value_hi-m_value, m_value-value_lo);
354   -
355   - // Write errors
356   - if (logTerse()) {
357   - log << gammalib::parformat("Error from curvature");
358   - log << m_model_par->error();
359   - log << " " << m_model_par->unit() << std::endl;
360   - log << gammalib::parformat("Error from profile");
361   - log << std::abs(error*m_model_par->scale());
362   - log << " " << m_model_par->unit() << std::endl;
363   - log << gammalib::parformat("Negative profile error");
364   - log << std::abs(error_neg*m_model_par->scale());
365   - log << " " << m_model_par->unit() << std::endl;
366   - log << gammalib::parformat("Positive profile error");
367   - log << std::abs(error_pos*m_model_par->scale());
368   - log << " " << m_model_par->unit() << std::endl;
369   - }
  304 + double error = 0.5 * (value_hi - value_lo);
  305 + double error_neg = m_value - value_lo;
  306 + double error_pos = value_hi - m_value;
  307 + double error_value = std::abs(error*m_model_par->scale());
  308 + double error_value_neg = std::abs(error_neg*m_model_par->scale());
  309 + double error_value_pos = std::abs(error_pos*m_model_par->scale());
  310 +
  311 + // Write results into logger
  312 + std::string unit = " " + m_model_par->unit();
  313 + log_value(NORMAL, "Lower parameter factor", value_lo);
  314 + log_value(NORMAL, "Upper parameter factor", value_hi);
  315 + log_value(NORMAL, "Error from curvature",
  316 + gammalib::str(m_model_par->error()) + unit);
  317 + log_value(NORMAL, "Error from profile",
  318 + gammalib::str(error_value) + unit);
  319 + log_value(NORMAL, "Negative profile error",
  320 + gammalib::str(error_value_neg) + unit);
  321 + log_value(NORMAL, "Positive profile error",
  322 + gammalib::str(error_value_pos) + unit);
370 323  
371 324 // Save error result
372 325 model->at(i).factor_error(error);
... ... @@ -381,13 +334,8 @@ void cterror::run(void)
381 334 // Recover optimizer
382 335 m_opt = best_opt;
383 336  
384   - // Restore energy dispersion flag for all CTA observations
385   - for (int i = 0; i < m_obs.size(); ++i) {
386   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
387   - if (obs != NULL) {
388   - obs->response()->apply_edisp(save_edisp[i]);
389   - }
390   - }
  337 + // Restore energy dispersion flags of all CTA observations
  338 + restore_edisp(m_obs, save_edisp);
391 339  
392 340 // Return
393 341 return;
... ... @@ -402,27 +350,26 @@ void cterror::run(void)
402 350 void cterror::save(void)
403 351 {
404 352 // Write header
405   - if (logTerse()) {
406   - log << std::endl;
407   - log.header1("Save results");
408   - }
  353 + log_header1(TERSE, "Save results");
409 354  
410 355 // Get output filename
411 356 m_outmodel = (*this)["outmodel"].filename();
412 357  
413   - // Save only if filename is non-empty
  358 + // Save only if filename is not empty and not "NONE"
414 359 if (!m_outmodel.is_empty() &&
415 360 gammalib::toupper(m_outmodel.url()) != "NONE") {
416 361  
417 362 // Log filename
418   - if (logTerse()) {
419   - log << "Save errors into file \""+m_outmodel+"\".";
420   - log << std::endl;
421   - }
  363 + log_value(NORMAL, "Model definition file", m_outmodel.url());
422 364  
423 365 // Write results out as XML model
424   - m_obs.models().save(m_outmodel);
  366 + m_obs.models().save(m_outmodel.url());
  367 +
  368 + }
425 369  
  370 + // ... otherwise signal that file was not save
  371 + else {
  372 + log_value(NORMAL, "Model definition file", "NONE");
426 373 }
427 374  
428 375 // Return
... ... @@ -448,6 +395,7 @@ void cterror::init_members(void)
448 395 m_tol = 1.0e-3;
449 396 m_max_iter = 50;
450 397 m_value = 0.0;
  398 + m_chatter = static_cast<GChatter>(2);
451 399  
452 400 // Initialise protected members
453 401 m_obs.clear();
... ... @@ -480,6 +428,7 @@ void cterror::copy_members(const cterror&amp; app)
480 428 m_tol = app.m_tol;
481 429 m_max_iter = app.m_max_iter;
482 430 m_apply_edisp = app.m_apply_edisp;
  431 + m_chatter = app.m_chatter;
483 432  
484 433 // Copy protected members
485 434 m_obs = app.m_obs;
... ... @@ -559,6 +508,7 @@ void cterror::get_parameters(void)
559 508 // Read other parameters
560 509 m_tol = (*this)["tol"].real();
561 510 m_max_iter = (*this)["max_iter"].integer();
  511 + m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
562 512  
563 513 // Read ahead parameters that are only needed when the tool gets
564 514 // executed
... ... @@ -606,9 +556,7 @@ double cterror::error_bisection(const double&amp; min, const double&amp; max)
606 556 "consider setting the minimum parameter "
607 557 "value to a lower value, and re-run "
608 558 "cterror.";
609   - if (logTerse()) {
610   - log << msg;
611   - }
  559 + log_string(TERSE, msg);
612 560 break;
613 561 }
614 562 else if (m_model_par->factor_max() - wrk_max < m_tol) {
... ... @@ -618,9 +566,7 @@ double cterror::error_bisection(const double&amp; min, const double&amp; max)
618 566 "consider setting the maximum parameter "
619 567 "value to a higher value, and re-run "
620 568 "cterror.";
621   - if (logTerse()) {
622   - log << msg;
623   - }
  569 + log_string(TERSE, msg);
624 570 break;
625 571 }
626 572 else {
... ... @@ -639,17 +585,9 @@ double cterror::error_bisection(const double&amp; min, const double&amp; max)
639 585 // Calculate function value
640 586 double eval_mid = evaluate(mid);
641 587  
642   - // Log interval
643   - if (logExplicit()) {
644   - log << gammalib::parformat(" Iteration "+gammalib::str(iter));
645   - log << "[";
646   - log << wrk_min;
647   - log << ", ";
648   - log << wrk_max;
649   - log << "] ";
650   - log << eval_mid;
651   - log << std::endl;
652   - }
  588 + // Write interval into logger
  589 + log_value(EXPLICIT, " Iteration "+gammalib::str(iter),
  590 + "["+gammalib::str(wrk_min)+", "+gammalib::str(wrk_max)+"]");
653 591  
654 592 // Check for convergence inside tolerance
655 593 if (std::abs(eval_mid) < m_tol) {
... ...
src/cterror/cterror.hpp
  View file @115eb8a
... ... @@ -35,7 +35,7 @@
35 35  
36 36 /* __Definitions _________________________________________________________ */
37 37 #define CTERROR_NAME "cterror"
38   -#define CTERROR_VERSION "1.1.0"
  38 +#define CTERROR_VERSION "1.2.0"
39 39  
40 40  
41 41 /***********************************************************************//**
... ... @@ -80,6 +80,7 @@ protected:
80 80 int m_max_iter; //!< Maximum number of iterations
81 81 double m_value; //!< Parameter value
82 82 bool m_apply_edisp; //!< Apply energy dispersion?
  83 + GChatter m_chatter; //!< Chattiness
83 84  
84 85 // Protected members
85 86 GObservations m_obs; //!< Observation container
... ... @@ -87,7 +88,6 @@ protected:
87 88 GModelPar* m_model_par; //!< Pointer to model parameter
88 89 double m_best_logL; //!< Best fit log likelihood of given model
89 90 GOptimizerLM m_opt; //!< Optimizer
90   -
91 91 };
92 92  
93 93  
... ...
src/ctlike/ctlike.cpp
  View file @115eb8a
... ... @@ -215,24 +215,12 @@ void ctlike::run(void)
215 215 // Get 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   - log.header1(gammalib::number("Observation",m_obs.size()));
234   - log << m_obs.print(m_chatter) << std::endl;
235   - }
  222 + // Write input observation container into logger
  223 + log_observations(NORMAL, m_obs, "Input observation");
236 224  
237 225 // Optimize model parameters using LM optimizer
238 226 optimize_lm();
... ... @@ -303,29 +291,16 @@ void ctlike::run(void)
303 291 }
304 292  
305 293 // Write results into logger
306   - if (logTerse()) {
307   - log << std::endl;
308   - log.header1("Maximum likelihood optimisation results");
309   - log << m_opt->print(m_chatter) << std::endl;
310   - log << gammalib::parformat("Maximum log likelihood");
311   - log << gammalib::str(m_logL,3) << std::endl;
312   - log << gammalib::parformat("Observed events (Nobs)");
313   - log << gammalib::str(num_events, 3) << std::endl;
314   - log << gammalib::parformat("Predicted events (Npred)");
315   - log << gammalib::str(npred, 3);
316   - log << " (Nobs - Npred = ";
317   - log << gammalib::str(num_events-npred);
318   - log << ")" << std::endl;
319   - log << m_obs.models().print(m_chatter) << std::endl;
320   - }
  294 + log_header1(NORMAL, "Maximum likelihood optimisation results");
  295 + log_string(NORMAL, m_opt->print(m_chatter));
  296 + log_value(NORMAL, "Maximum log likelihood", gammalib::str(m_logL,3));
  297 + log_value(NORMAL, "Observed events (Nobs)", gammalib::str(num_events,3));
  298 + log_value(NORMAL, "Predicted events (Npred)", gammalib::str(npred,3)+
  299 + " (Nobs - Npred = "+gammalib::str(num_events-npred)+")");
  300 + log_string(NORMAL, m_obs.models().print(m_chatter));
321 301  
322   - // Restore energy dispersion flag for all CTA observations
323   - for (int i = 0; i < m_obs.size(); ++i) {
324   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
325   - if (obs != NULL) {
326   - obs->response()->apply_edisp(save_edisp[i]);
327   - }
328   - }
  302 + // Restore energy dispersion flags of all CTA observations
  303 + restore_edisp(m_obs, save_edisp);
329 304  
330 305 // Return
331 306 return;
... ... @@ -340,10 +315,7 @@ void ctlike::run(void)
340 315 void ctlike::save(void)
341 316 {
342 317 // Write header
343   - if (logTerse()) {
344   - log << std::endl;
345   - log.header1("Save results");
346   - }
  318 + log_header1(TERSE, "Save results");
347 319  
348 320 // Get output filename
349 321 m_outmodel = (*this)["outmodel"].filename();
... ... @@ -353,16 +325,18 @@ void ctlike::save(void)
353 325 gammalib::toupper(m_outmodel.url()) != "NONE") {
354 326  
355 327 // Log filename
356   - if (logTerse()) {
357   - log << "Save fitted parameters into file \""+m_outmodel+"\".";
358   - log << std::endl;
359   - }
  328 + log_value(NORMAL, "Model definition file", m_outmodel.url());
360 329  
361 330 // Write results out as XML model
362 331 m_obs.models().save(m_outmodel);
363 332  
364 333 }
365 334  
  335 + // ... otherwise signal that file was not save
  336 + else {
  337 + log_value(NORMAL, "Model definition file", "NONE");
  338 + }
  339 +
366 340 // Return
367 341 return;
368 342 }
... ... @@ -434,7 +408,7 @@ void ctlike::get_parameters(void)
434 408 }
435 409  
436 410 // Set optimizer logger
437   - if (logTerse()) {
  411 + if (logNormal()) {
438 412 static_cast<GOptimizerLM*>(m_opt)->logger(&log);
439 413 }
440 414 else {
... ... @@ -456,12 +430,9 @@ void ctlike::get_parameters(void)
456 430 ***************************************************************************/
457 431 void ctlike::optimize_lm(void)
458 432 {
459   - // Write Header for optimization and indent for optimizer logging
460   - if (logTerse()) {
461   - log << std::endl;
462   - log.header1("Maximum likelihood optimisation");
463   - log.indent(1);
464   - }
  433 + // Write header
  434 + log_header1(TERSE, "Maximum likelihood optimisation");
  435 + log.indent(1);
465 436  
466 437 // Compute number of fitted parameters
467 438 int nfit = 0;
... ... @@ -476,14 +447,9 @@ void ctlike::optimize_lm(void)
476 447  
477 448 // Notify if all parameters are fixed
478 449 if (nfit == 0) {
479   - if (logTerse()) {
480   - log << "WARNING: All model parameters are fixed!";
481   - log << std::endl;
482   - log << " ctlike will proceed without fitting parameters.";
483   - log << std::endl;
484   - log << " All curvature matrix elements will be zero.";
485   - log << std::endl;
486   - }
  450 + log_string(TERSE, "WARNING: All model parameters are fixed!");
  451 + log_string(TERSE, " ctlike will proceed without fitting parameters.");
  452 + log_string(TERSE, " All curvature matrix elements will be zero.");
487 453 }
488 454  
489 455 // Perform LM optimization
... ... @@ -494,11 +460,8 @@ void ctlike::optimize_lm(void)
494 460  
495 461 // Dump new header
496 462 log.indent(0);
497   - if (logTerse()) {
498   - log << std::endl;
499   - log.header1("Maximum likelihood re-optimisation");
500   - log.indent(1);
501   - }
  463 + log_header1(TERSE, "Maximum likelihood re-optimisation");
  464 + log.indent(1);
502 465  
503 466 // Optimise again
504 467 m_obs.optimize(*m_opt);
... ... @@ -506,13 +469,10 @@ void ctlike::optimize_lm(void)
506 469 }
507 470  
508 471 // Optionally show curvature matrix
509   - if (logNormal()) {
510   - log << std::endl;
511   - log.header1("Curvature matrix");
512   - log.indent(1);
513   - log << *(const_cast<GObservations::likelihood&>(m_obs.function()).curvature());
514   - log << std::endl;
515   - }
  472 + log_header1(EXPLICIT, "Curvature matrix");
  473 + log.indent(1);
  474 + log_string(EXPLICIT, (const_cast<GObservations::likelihood&>
  475 + (m_obs.function()).curvature())->print());
516 476  
517 477 // Compute errors
518 478 m_obs.errors(*m_opt);
... ... @@ -537,11 +497,8 @@ void ctlike::optimize_lm(void)
537 497 double ctlike::reoptimize_lm(void)
538 498 {
539 499 // Write Header for optimization and indent for optimizer logging
540   - if (logTerse()) {
541   - log << std::endl;
542   - log.header1("Maximum likelihood re-optimisation");
543   - log.indent(1);
544   - }
  500 + log_header1(TERSE, "Maximum likelihood re-optimisation");
  501 + log.indent(1);
545 502  
546 503 // Create a clone of the optimizer for the re-optimisation
547 504 GOptimizer* opt = m_opt->clone();
... ... @@ -559,11 +516,8 @@ double ctlike::reoptimize_lm(void)
559 516  
560 517 // Write optimization results
561 518 log.indent(0);
562   - if (logTerse()) {
563   - log << std::endl;
564   - log.header1("Maximum likelihood re-optimisation results");
565   - log << *opt << std::endl;
566   - }
  519 + log_header1(NORMAL, "Maximum likelihood re-optimisation results");
  520 + log_string(NORMAL, opt->print(m_chatter));
567 521  
568 522 // Return
569 523 return (logL);
... ...
src/support/ctool.cpp
  View file @115eb8a
... ... @@ -52,6 +52,8 @@
52 52 #define G_SETUP_OBSERVATION "ctool::setup_observations(GObservations&)"
53 53 #define G_GET_MEAN_POINTING "ctool::get_mean_pointing(GObservations&)"
54 54 #define G_CREATE_EBOUNDS "ctool::create_ebounds()"
  55 +#define G_RESTORE_EDISP "ctool::restore_edisp(GObservations&, "\
  56 + "std::vector<bool>&)"
55 57 #define G_SET_OBS_RESPONSE "ctool::set_obs_response(GCTAObservation*)"
56 58 #define G_PROVIDE_HELP "ctool::provide_help()"
57 59  
... ... @@ -994,7 +996,7 @@ void ctool::set_response(GObservations&amp; obs)
994 996 // Is this observation a CTA observation?
995 997 GCTAObservation* cta = dynamic_cast<GCTAObservation*>(obs[i]);
996 998  
997   - // Yes ...
  999 + // Yes, the set the response if the observation does not yet have one
998 1000 if (cta != NULL) {
999 1001  
1000 1002 // Set response if we don't have one
... ... @@ -1012,6 +1014,87 @@ void ctool::set_response(GObservations&amp; obs)
1012 1014  
1013 1015  
1014 1016 /***********************************************************************//**
  1017 + * @brief Set energy dispersion to CTA observations
  1018 + *
  1019 + * @param[in] obs Observation container.
  1020 + * @param[in] edisp Requested energy dispersion flag.
  1021 + * @return Vector of old energy dispersion flags.
  1022 + *
  1023 + * Applies energy dispersion to all CTA observations. The method returns a
  1024 + * vector with the old energy dispersion flags.
  1025 + ***************************************************************************/
  1026 +std::vector<bool> ctool::set_edisp(GObservations& obs, const bool& edisp) const
  1027 +{
  1028 + // Allocate vector for energy dispersion flags of all observation in
  1029 + // the container
  1030 + std::vector<bool> old_edisp(obs.size(), false);
  1031 +
  1032 + // Loop over all observations in observation container
  1033 + for (int i = 0; i < obs.size(); ++i) {
  1034 +
  1035 + // Is this observation a CTA observation?
  1036 + GCTAObservation* cta = dynamic_cast<GCTAObservation*>(obs[i]);
  1037 +
  1038 + // Yes, then set the energy dispersion flag
  1039 + if (cta != NULL) {
  1040 +
  1041 + // Save old energy dispersion flag
  1042 + old_edisp[i] = cta->response()->apply_edisp();
  1043 +
  1044 + // Set energy dispersion flag according to input parameter
  1045 + cta->response()->apply_edisp(edisp);
  1046 +
  1047 + } // endif: observation was CTA observation
  1048 +
  1049 + } // endfor: loop over all observations in container
  1050 +
  1051 + // Return old energy dispersion flags
  1052 + return old_edisp;
  1053 +}
  1054 +
  1055 +
  1056 +/***********************************************************************//**
  1057 + * @brief Restore energy dispersion flags of CTA observations
  1058 + *
  1059 + * @param[in] obs Observation container.
  1060 + * @param[in] edisp Vector of energy dispersion flags.
  1061 + *
  1062 + * Restores energy dispersion flags that have previously been returned by
  1063 + * the set_edisp() method. The number of observations between the calls
  1064 + * of the set_edisp() and this method has to be the same.
  1065 + ***************************************************************************/
  1066 +void ctool::restore_edisp(GObservations& obs, const std::vector<bool>& edisp) const
  1067 +{
  1068 + // Check that energy dispersion flag vector has the same size as the
  1069 + // observation container
  1070 + if (edisp.size() != obs.size()) {
  1071 + std::string msg = "Size of energy dispersion flag vector ("+
  1072 + gammalib::str(edisp.size())+" elements) is "
  1073 + "incompatible with number of observations in the "
  1074 + "observation container ("+gammalib::str(obs.size())+
  1075 + " observation).";
  1076 + throw GException::invalid_value(G_RESTORE_EDISP, msg);
  1077 + }
  1078 +
  1079 + // Loop over all observations in observation container
  1080 + for (int i = 0; i < obs.size(); ++i) {
  1081 +
  1082 + // Is this observation a CTA observation?
  1083 + GCTAObservation* cta = dynamic_cast<GCTAObservation*>(obs[i]);
  1084 +
  1085 + // Yes, then restore the energy dispersion flag
  1086 + if (cta != NULL) {
  1087 + cta->response()->apply_edisp(edisp[i]);
  1088 + }
  1089 +
  1090 + } // endfor: loop over all observations in container
  1091 +
  1092 + // Return
  1093 + return;
  1094 +}
  1095 +
  1096 +
  1097 +/***********************************************************************//**
1015 1098 * @brief Set response for CTA observation
1016 1099 *
1017 1100 * @param[in,out] obs CTA observation
... ...
src/support/ctool.hpp
  View file @115eb8a
... ... @@ -94,14 +94,17 @@ protected:
94 94 const std::string& what = "Model");
95 95  
96 96 // Protected support methods
97   - void set_response(GObservations& obs);
98   - void set_obs_response(GCTAObservation* obs);
99   - void set_obs_bounds(GObservations& obs);
100   - GSkyDir get_mean_pointing(const GObservations& obs);
101   - size_t get_current_rss(void);
102   - std::string get_obs_header(const GObservation* obs);
103   - GEnergies insert_energy_boundaries(const GEnergies& energies,
104   - const GCTAObservation& obs);
  97 + void set_response(GObservations& obs);
  98 + std::vector<bool> set_edisp(GObservations& obs, const bool& edisp) const;
  99 + void restore_edisp(GObservations& obs,
  100 + const std::vector<bool>& edisp) const;
  101 + void set_obs_response(GCTAObservation* obs);
  102 + void set_obs_bounds(GObservations& obs);
  103 + GSkyDir get_mean_pointing(const GObservations& obs);
  104 + size_t get_current_rss(void);
  105 + std::string get_obs_header(const GObservation* obs);
  106 + GEnergies insert_energy_boundaries(const GEnergies& energies,
  107 + const GCTAObservation& obs);
105 108  
106 109 // Protected warning strings
107 110 std::string warn_too_few_energies(const GEnergies& energies) const;
... ...
test/test_cterror.py
  View file @115eb8a
... ... @@ -101,14 +101,31 @@ class Test(test):
101 101 """
102 102 Test cterror from Python
103 103 """
104   - # Set-up cterror
  104 + # Allocate cterror
105 105 error = ctools.cterror()
  106 +
  107 + # Check that empty cterror tool holds an empty observation and
  108 + # optimizer
  109 + self._check_obs(error.obs(), nobs=0, nmodels=0)
  110 + self._check_opt(error.opt())
  111 +
  112 + # Check that saving saves an empty model definition file
  113 + error['outmodel'] = 'cterror_py0.xml'
  114 + error['logfile'] = 'cterror_py0.log'
  115 + error.logFileOpen()
  116 + error.save()
  117 + self._check_result_file('cterror_py0.xml', nmodels=0)
  118 +
  119 + # Check that clearing does not lead to an exception or segfault
  120 + error.clear()
  121 +
  122 + # Now set cterror parameters
106 123 error['inobs'] = self._events
107 124 error['inmodel'] = self._model
108   - error['outmodel'] = 'cterror_py1.xml'
109 125 error['srcname'] = 'Crab'
110 126 error['caldb'] = self._caldb
111 127 error['irf'] = self._irf
  128 + error['outmodel'] = 'cterror_py1.xml'
112 129 error['logfile'] = 'cterror_py1.log'
113 130 error['chatter'] = 2
114 131  
... ... @@ -120,38 +137,130 @@ class Test(test):
120 137 # Check result file
121 138 self._check_result_file('cterror_py1.xml')
122 139  
  140 + # Copy cterror tool
  141 + cpy_error = error.copy()
  142 +
  143 + # Check observation and optimizer of copy
  144 + self._check_obs(cpy_error.obs())
  145 + self._check_opt(cpy_error.opt())
  146 +
  147 + # Execute copy of cterror tool again, now with a higher chatter
  148 + # level than before
  149 + cpy_error['outmodel'] = 'cterror_py2.xml'
  150 + cpy_error['logfile'] = 'cterror_py2.log'
  151 + cpy_error['chatter'] = 3
  152 + cpy_error.logFileOpen() # Needed to get a new log file
  153 + cpy_error.execute()
  154 +
  155 + # Check result file
  156 + self._check_result_file('cterror_py2.xml')
  157 +
  158 + # Now clear copy of cterror tool
  159 + cpy_error.clear()
  160 +
  161 + # Check that the cleared copy has also cleared the observations
  162 + self._check_obs(cpy_error.obs(), nobs=0, nmodels=0)
  163 +
123 164 # Return
124 165 return
125 166  
126 167 # Check result file
127   - def _check_result_file(self, filename):
  168 + def _check_result_file(self, filename, nmodels=2):
128 169 """
129 170 Check result file
  171 +
  172 + Parameters
  173 + ----------
  174 + filename : str
  175 + Model definition XML file
  176 + nmodels : int, optional
  177 + Expected number of models
  178 + """
  179 + # Open models from result file
  180 + models = gammalib.GModels(filename)
  181 +
  182 + # Check models
  183 + self._check_models(models, nmodels=nmodels)
  184 +
  185 + # Return
  186 + return
  187 +
  188 + # Check observation and models
  189 + def _check_obs(self, obs, nobs=1, nmodels=2):
  190 + """
  191 + Check observation and models
  192 +
  193 + Parameters
  194 + ----------
  195 + obs : `~gammalib.GObservations`
  196 + Models
  197 + nobs : int, optional
  198 + Expected number of observations
  199 + nmodels : int, optional
  200 + Expected number of models
  201 + """
  202 + # Check number of observations
  203 + self.test_value(obs.size(), nobs, 'Check number of observations')
  204 +
  205 + # Check models
  206 + self._check_models(obs.models(), nmodels=nmodels)
  207 +
  208 + # Return
  209 + return
  210 +
  211 + # Check models
  212 + def _check_models(self, models, nmodels=2):
  213 + """
  214 + Check observation and models
  215 +
  216 + Parameters
  217 + ----------
  218 + models : `~gammalib.GModels`
  219 + Models
  220 + nmodels : int, optional
  221 + Expected number of models
  222 + """
  223 + # Check number of models
  224 + self.test_value(models.size(), nmodels, 'Check number of models')
  225 +
  226 + # If there are models in the container then check them
  227 + if nmodels > 0:
  228 + self.test_value(models['Crab']['Prefactor'].value(), 1.58907e-16,
  229 + 1.0e-3, 'Check fitted Crab Prefactor')
  230 + self.test_value(models['Crab']['Prefactor'].error(), 0.0529105e-16,
  231 + 1.0e-3, 'Check Crab Prefactor error')
  232 + self.test_value(models['Crab']['Index'].value(), -2.43549,
  233 + 1.0e-3, 'Check fitted Crab Index')
  234 + self.test_value(models['Crab']['Index'].error(), 0.027804,
  235 + 1.0e-3, 'Check Crab Index error')
  236 + self.test_value(models['Background']['Prefactor'].value(), 61.6919e-6,
  237 + 1.0e-3, 'Check fitted background Prefactor')
  238 + self.test_value(models['Background']['Prefactor'].error(), 1.49438e-6,
  239 + 1.0e-3, 'Check background Prefactor error')
  240 + self.test_value(models['Background']['Index'].value(), -2.20535,
  241 + 1.0e-3, 'Check fitted background Index')
  242 + self.test_value(models['Background']['Index'].error(), 0.0113269,
  243 + 1.0e-3, 'Check background Index error')
  244 + self.test_value(models['Background']['Sigma'].value(), 3.04252,
  245 + 1.0e-3, 'Check fitted background Sigma')
  246 + self.test_value(models['Background']['Sigma'].error(), 0.0307008,
  247 + 1.0e-3, 'Check background Sigma error')
  248 +
  249 + # Return
  250 + return
  251 +
  252 + # Check optimizer
  253 + def _check_opt(self, opt):
  254 + """
  255 + Check optimizer
  256 +
  257 + Parameters
  258 + ----------
  259 + opt : `~gammalib.GOptimizerLM`
  260 + Optimizer
130 261 """
131   - # Open result file
132   - result = gammalib.GModels(filename)
133   -
134   - # Check results
135   - self.test_value(result['Crab']['Prefactor'].value(), 1.58907e-16,
136   - 1.0e-3, 'Check fitted Crab Prefactor')
137   - self.test_value(result['Crab']['Prefactor'].error(), 0.0529105e-16,
138   - 1.0e-3, 'Check Crab Prefactor error')
139   - self.test_value(result['Crab']['Index'].value(), -2.43549,
140   - 1.0e-3, 'Check fitted Crab Index')
141   - self.test_value(result['Crab']['Index'].error(), 0.027804,
142   - 1.0e-3, 'Check Crab Index error')
143   - self.test_value(result['Background']['Prefactor'].value(), 61.6919e-6,
144   - 1.0e-3, 'Check fitted background Prefactor')
145   - self.test_value(result['Background']['Prefactor'].error(), 1.49438e-6,
146   - 1.0e-3, 'Check background Prefactor error')
147   - self.test_value(result['Background']['Index'].value(), -2.20535,
148   - 1.0e-3, 'Check fitted background Index')
149   - self.test_value(result['Background']['Index'].error(), 0.0113269,
150   - 1.0e-3, 'Check background Index error')
151   - self.test_value(result['Background']['Sigma'].value(), 3.04252,
152   - 1.0e-3, 'Check fitted background Sigma')
153   - self.test_value(result['Background']['Sigma'].error(), 0.0307008,
154   - 1.0e-3, 'Check background Sigma error')
  262 + # Check optimizer
  263 + self.test_value(opt.status(), 0, 'Check optimizer status')
155 264  
156 265 # Return
157 266 return
... ...
test/test_ctlike.py
  View file @115eb8a
... ... @@ -101,8 +101,25 @@ class Test(test):
101 101 """
102 102 Test ctlike from Python
103 103 """
104   - # Set-up ctlike
  104 + # Allocate ctlike
105 105 like = ctools.ctlike()
  106 +
  107 + # Check that empty ctlike tool holds an empty observation and
  108 + # optimizer
  109 + self._check_obs(like.obs(), nobs=0, nmodels=0)
  110 + self._check_opt(like.opt())
  111 +
  112 + # Check that saving saves an empty model definition file
  113 + like['outmodel'] = 'ctlike_py0.xml'
  114 + like['logfile'] = 'ctlike_py0.log'
  115 + like.logFileOpen()
  116 + like.save()
  117 + self._check_result_file('ctlike_py0.xml', nmodels=0)
  118 +
  119 + # Check that clearing does not lead to an exception or segfault
  120 + like.clear()
  121 +
  122 + # Now set ctlike parameters
106 123 like['inobs'] = self._events
107 124 like['inmodel'] = self._model
108 125 like['caldb'] = self._caldb
... ... @@ -119,46 +136,130 @@ class Test(test):
119 136 # Check result file
120 137 self._check_result_file('ctlike_py1.xml')
121 138  
122   - # Retrieve observation container, save and check model file
123   - obs = like.obs()
124   - obs.models().save('ctlike_py2.xml')
  139 + # Copy ctlike tool
  140 + cpy_like = like.copy()
  141 +
  142 + # Check observation and optimizer of copy
  143 + self._check_obs(cpy_like.obs())
  144 + self._check_opt(cpy_like.opt())
  145 +
  146 + # Execute copy of ctlike tool again, now with a higher chatter
  147 + # level than before
  148 + cpy_like['outmodel'] = 'ctlike_py2.xml'
  149 + cpy_like['logfile'] = 'ctlike_py2.log'
  150 + cpy_like['chatter'] = 3
  151 + cpy_like.logFileOpen() # Needed to get a new log file
  152 + cpy_like.execute()
  153 +
  154 + # Check result file
125 155 self._check_result_file('ctlike_py2.xml')
126 156  
127   - # Retrieve optimizer
128   - opt = like.opt()
  157 + # Now clear copy of ctlike tool
  158 + cpy_like.clear()
  159 +
  160 + # Check that the cleared copy has also cleared the observations
  161 + self._check_obs(cpy_like.obs(), nobs=0, nmodels=0)
129 162  
130 163 # Return
131 164 return
132 165  
133 166 # Check result file
134   - def _check_result_file(self, filename):
  167 + def _check_result_file(self, filename, nmodels=2):
135 168 """
136 169 Check result file
  170 +
  171 + Parameters
  172 + ----------
  173 + filename : str
  174 + Model definition XML file
  175 + nmodels : int, optional
  176 + Expected number of models
  177 + """
  178 + # Open models from result file
  179 + models = gammalib.GModels(filename)
  180 +
  181 + # Check models
  182 + self._check_models(models, nmodels=nmodels)
  183 +
  184 + # Return
  185 + return
  186 +
  187 + # Check observation and models
  188 + def _check_obs(self, obs, nobs=1, nmodels=2):
  189 + """
  190 + Check observation and models
  191 +
  192 + Parameters
  193 + ----------
  194 + obs : `~gammalib.GObservations`
  195 + Models
  196 + nobs : int, optional
  197 + Expected number of observations
  198 + nmodels : int, optional
  199 + Expected number of models
  200 + """
  201 + # Check number of observations
  202 + self.test_value(obs.size(), nobs, 'Check number of observations')
  203 +
  204 + # Check models
  205 + self._check_models(obs.models(), nmodels=nmodels)
  206 +
  207 + # Return
  208 + return
  209 +
  210 + # Check models
  211 + def _check_models(self, models, nmodels=2):
  212 + """
  213 + Check observation and models
  214 +
  215 + Parameters
  216 + ----------
  217 + models : `~gammalib.GModels`
  218 + Models
  219 + nmodels : int, optional
  220 + Expected number of models
  221 + """
  222 + # Check number of models
  223 + self.test_value(models.size(), nmodels, 'Check number of models')
  224 +
  225 + # If there are models in the container then check them
  226 + if nmodels > 0:
  227 + self.test_value(models['Crab']['Prefactor'].value(), 1.58907e-16,
  228 + 1.0e-3, 'Check fitted Crab Prefactor')
  229 + self.test_value(models['Crab']['Prefactor'].error(), 0.0526982e-16,
  230 + 1.0e-3, 'Check Crab Prefactor error')
  231 + self.test_value(models['Crab']['Index'].value(), -2.43549,
  232 + 1.0e-3, 'Check fitted Crab Index')
  233 + self.test_value(models['Crab']['Index'].error(), 0.0248116,
  234 + 1.0e-3, 'Check Crab Index error')
  235 + self.test_value(models['Background']['Prefactor'].value(), 61.6919e-6,
  236 + 1.0e-3, 'Check fitted background Prefactor')
  237 + self.test_value(models['Background']['Prefactor'].error(), 1.49438e-6,
  238 + 1.0e-3, 'Check background Prefactor error')
  239 + self.test_value(models['Background']['Index'].value(), -2.20535,
  240 + 1.0e-3, 'Check fitted background Index')
  241 + self.test_value(models['Background']['Index'].error(), 0.0113269,
  242 + 1.0e-3, 'Check background Index error')
  243 + self.test_value(models['Background']['Sigma'].value(), 3.04252,
  244 + 1.0e-3, 'Check fitted background Sigma')
  245 + self.test_value(models['Background']['Sigma'].error(), 0.0307008,
  246 + 1.0e-3, 'Check background Sigma error')
  247 +
  248 + # Return
  249 + return
  250 +
  251 + # Check optimizer
  252 + def _check_opt(self, opt):
  253 + """
  254 + Check optimizer
  255 +
  256 + Parameters
  257 + ----------
  258 + opt : `~gammalib.GOptimizerLM`
  259 + Optimizer
137 260 """
138   - # Open result file
139   - result = gammalib.GModels(filename)
140   -
141   - # Check results
142   - self.test_value(result['Crab']['Prefactor'].value(), 1.58907e-16,
143   - 1.0e-3, 'Check fitted Crab Prefactor')
144   - self.test_value(result['Crab']['Prefactor'].error(), 0.0526982e-16,
145   - 1.0e-3, 'Check Crab Prefactor error')
146   - self.test_value(result['Crab']['Index'].value(), -2.43549,
147   - 1.0e-3, 'Check fitted Crab Index')
148   - self.test_value(result['Crab']['Index'].error(), 0.0248116,
149   - 1.0e-3, 'Check Crab Index error')
150   - self.test_value(result['Background']['Prefactor'].value(), 61.6919e-6,
151   - 1.0e-3, 'Check fitted background Prefactor')
152   - self.test_value(result['Background']['Prefactor'].error(), 1.49438e-6,
153   - 1.0e-3, 'Check background Prefactor error')
154   - self.test_value(result['Background']['Index'].value(), -2.20535,
155   - 1.0e-3, 'Check fitted background Index')
156   - self.test_value(result['Background']['Index'].error(), 0.0113269,
157   - 1.0e-3, 'Check background Index error')
158   - self.test_value(result['Background']['Sigma'].value(), 3.04252,
159   - 1.0e-3, 'Check fitted background Sigma')
160   - self.test_value(result['Background']['Sigma'].error(), 0.0307008,
161   - 1.0e-3, 'Check background Sigma error')
  261 + # Check optimizer
  262 + self.test_value(opt.status(), 0, 'Check optimizer status')
162 263  
163 264 # Return
164 265 return
... ...