Francois Brun / ctools

Commit 411b731669f58bb0b28c6b1874d3cb89df2b2311

Authored by Jürgen Knödlseder
1 parent 24cd52b8

Integrate ctulimit tool.

Inline Side-by-side
Showing 7 changed files with 346 additions and 284 deletions
doc/source/reference_manual/ctulimit.rst
  View file @411b731
... ... @@ -3,14 +3,15 @@
3 3 ctulimit
4 4 ===========
5 5  
6   -Computes upper limit.
  6 +Computes upper limit for a specific sky model component.
7 7  
8 8  
9 9 Synopsis
10 10 --------
11 11  
12   -Calculates the upper limit flux of a given source and stores the value in ascii file.
13   -The input model has to be optimised for ctulimit to converge.
  12 +Computes the upper limit flux of a given source and stores the value
  13 +in an ASCII file.
  14 +
14 15  
15 16 General parameters
16 17 ------------------
... ... @@ -39,26 +40,29 @@ General parameters
39 40 ``outfile [file]``
40 41 Output ASCII file name.
41 42  
42   -``(cl = 0.95) [real]``
43   - Confidence Level of upper limit
  43 +``(confidence = 0.95) [real]``
  44 + Confidence Level of upper limit.
44 45  
45 46 ``(sigma_min = 0.0) [real]``
46   - Minimum boundary to start searching for upper limit value. Number of standard deviations above best fit value
  47 + Minimum boundary to start searching for upper limit value.
  48 + Number of standard deviations above best fit value
47 49  
48 50 ``(sigma_max = 10.0) [real]``
49   - Maximum boundary to start searching for upper limit value. Number of standard deviations above best fit value
  51 + Maximum boundary to start searching for upper limit value.
  52 + Number of standard deviations above best fit value
50 53  
51 54 ``(emin = 1.0) [real])``
52   - Minimum energy of flux upper limit (in TeV).
  55 + Minimum energy of flux limits (in TeV).
53 56  
54 57 ``(emax = 100) [real]``
55   - Maximum energy of flux upper limit (in TeV).
  58 + Maximum energy of flux limits (in TeV).
56 59  
57 60 ``(tol = 1e-5) [real]``
58   - Computation tolerance for minimum
  61 + Computation tolerance.
59 62  
60 63 ``(max_iter = 50) [integer]``
61   - Maximum number of iterations before throwing an exception
  64 + Maximum number of iterations before stopping the upper
  65 + limit computations.
62 66  
63 67  
64 68 Standard parameters
... ... @@ -85,7 +89,7 @@ Standard parameters
85 89 ``(mode = ql) [string]``
86 90 Mode of automatic parameters (default is "ql", i.e. "query and learn").
87 91  
88   -``(logfile = ctskymap.log) [string]``
  92 +``(logfile = ctulimit.log) [string]``
89 93 Name of log file.
90 94  
91 95  
... ...
pyext/ctulimit.i
  View file @411b731
1 1 /***************************************************************************
2   - * ctulimit - ctulimit calculation tool *
  2 + * ctulimit - Upper limit calculation tool *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2015 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file ctulimit
23   - * @brief upper limit calculation tool
  23 + * @brief Upper limit calculation tool
24 24 * @author Michael Mayer
25 25 */
26 26 %{
... ... @@ -32,7 +32,7 @@
32 32 /***********************************************************************//**
33 33 * @class ctulimit
34 34 *
35   - * @brief upper limit calculation tool
  35 + * @brief Upper limit calculation tool
36 36 ***************************************************************************/
37 37 class ctulimit : public ctool {
38 38 public:
... ... @@ -48,14 +48,14 @@ public:
48 48 void run(void);
49 49 void save(void);
50 50 const GObservations& obs(void) const;
51   - const double& diff_ulimit(void) const;
52   - const double& flux_ulimit(void) const;
53   - const double& eflux_ulimit(void) const;
  51 + const double& diff_ulimit(void) const;
  52 + const double& flux_ulimit(void) const;
  53 + const double& eflux_ulimit(void) const;
54 54 };
55 55  
56 56  
57 57 /***********************************************************************//**
58   - * @brief upper limit calculation tool Python extensions
  58 + * @brief Upper limit calculation tool Python extensions
59 59 ***************************************************************************/
60 60 %extend ctulimit {
61 61 ctulimit copy() {
... ...
src/ctulimit/ctulimit.cpp
  View file @411b731
1 1 /***************************************************************************
2   - * ctulimit - upper limit calculation tool *
  2 + * ctulimit - Upper limit calculation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Michael Mayer *
  4 + * copyright (C) 2015 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file ctulimit.hpp
23   - * @brief upper limit calculation tool interface implementation
  23 + * @brief Upper limit calculation tool interface implementation
24 24 * @author Michael Mayer
25 25 */
26 26  
... ... @@ -34,9 +34,11 @@
34 34 #include "GOptimizer.hpp"
35 35  
36 36 /* __ Method name definitions ____________________________________________ */
37   -#define G_GET_PARAMETERS "ctulimit::get_parameters()"
38   -#define G_UL_BISECTION "ctulimit::ul_bisection(const double&, const double&)"
39   -#define G_EVALUATE "ctulimit::evaluate(const double&)"
  37 +#define G_GET_PARAMETERS "ctulimit::get_parameters()"
  38 +#define G_GET_MODEL_PARAMETER "ctulimit::get_model_parameter()"
  39 +#define G_UL_BISECTION "ctulimit::ul_bisection(double&, double&)"
  40 +#define G_EVALUATE "ctulimit::evaluate(double&)"
  41 +
40 42 /* __ Debug definitions __________________________________________________ */
41 43  
42 44 /* __ Coding definitions _________________________________________________ */
... ... @@ -70,7 +72,7 @@ ctulimit::ctulimit(void) : ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
70 72 * observations container.
71 73 ***************************************************************************/
72 74 ctulimit::ctulimit(const GObservations& obs) :
73   - ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
  75 + ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
74 76 {
75 77 // Initialise members
76 78 init_members();
... ... @@ -91,7 +93,7 @@ ctulimit::ctulimit(const GObservations& obs) :
91 93 * @param[in] argv Array of command line arguments.
92 94 ***************************************************************************/
93 95 ctulimit::ctulimit(int argc, char *argv[]) :
94   - ctool(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv)
  96 + ctool(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv)
95 97 {
96 98 // Initialise members
97 99 init_members();
... ... @@ -227,80 +229,97 @@ void ctulimit::run(void)
227 229 log << m_obs << std::endl;
228 230 }
229 231  
  232 + // Save original models
230 233 GModels models_orig = m_obs.models();
231 234  
232   -
233   - m_bestloglike = m_obs.logL();
234   - if (m_bestloglike == 0.0) {
  235 + // Save original log-likelihood. If the value is zero it has never been
  236 + // computed hence we compute it now.
  237 + m_best_logL = m_obs.logL();
  238 + if (m_best_logL == 0.0) {
235 239  
236 240 // Write header
237 241 if (logTerse()) {
238 242 log << std::endl;
239   - log.header1("Recompute best-fit likelihood");
  243 + log.header1("Compute best-fit likelihood");
240 244 }
241 245  
242 246 // Reoptimize if likelihood was not given before
243 247 GOptimizerLM* opt = new GOptimizerLM();
244 248 m_obs.optimize(*opt);
245 249 m_obs.errors(*opt);
246   - m_bestloglike = m_obs.logL();
247   -
248   - } //endif: likelihood was 0.0
249   -
250   - // Store optimized models and pointer to skymodel of interest
251   - m_models = m_obs.models();
252   - m_skymodel = dynamic_cast<GModelSky*>(m_models[m_srcname])->clone();
253   -
254   - // Get reference energy for differential upper limit
255   - GEnergy eref = GEnergy(m_eref, "TeV");
  250 + m_best_logL = m_obs.logL();
256 251  
257   - // Create energy range for flux limits
258   - GEnergy emin = GEnergy(m_emin, "TeV");
259   - GEnergy emax = GEnergy(m_emax, "TeV");
  252 + // Write optimised model into logger
  253 + if (logTerse()) {
  254 + log << m_obs.models() << std::endl;
  255 + }
260 256  
261   - // Get value, scale and error of parameter
262   - double value = (*m_skymodel->spectral())[0].value();
263   - double error = (*m_skymodel->spectral())[0].error();
264   - double scale = (*m_skymodel->spectral())[0].scale();
  257 + } // endif: likelihood was zero
265 258  
266   - // Compute starting boundaries
267   - double parmin = value / scale + m_sigma_min * error / scale;
268   - double parmax = value / scale + m_sigma_max * error / scale;
  259 + // Compute parameter bracketing
  260 + double value = m_model_par->value();
  261 + double error = m_model_par->error();
  262 + double parmin = value + m_sigma_min * error;
  263 + double parmax = value + m_sigma_max * error;
269 264  
270 265 // Write header
271 266 if (logTerse()) {
272 267 log << std::endl;
273 268 log.header1("Compute upper limit");
274   - log << "Searching for upper limit between ";
275   - log << parmin * scale;
276   - log << " and ";
277   - log << parmax * scale;
278   - log << ": "<<std::endl;
  269 + log << gammalib::parformat("Model name");
  270 + log << m_skymodel->name() << std::endl;
  271 + log << gammalib::parformat("Parameter name");
  272 + log << m_model_par->name() << std::endl;
  273 + log << gammalib::parformat("Initial parameter range");
  274 + log << "[";
  275 + log << parmin;
  276 + log << ", ";
  277 + log << parmax;
  278 + log << "]" << std::endl;
279 279 }
280 280  
281   - // compute upper limit
282   - ulimit_bisection(parmin, parmax, scale);
  281 + // Compute upper limit
  282 + ulimit_bisection(parmin, parmax);
283 283  
284   - m_diff_ulimit = m_skymodel->spectral()->eval(eref, GTime());
285   - m_flux_ulimit = m_skymodel->spectral()->flux(emin, emax);
  284 + // Write final parameter range
  285 + if (logTerse()) {
  286 + log << gammalib::parformat("Final parameter");
  287 + log << m_model_par->value() << std::endl;
  288 + }
  289 +
  290 + // Get reference energy for differential upper limit
  291 + GEnergy eref = GEnergy(m_eref, "TeV");
  292 +
  293 + // Create energy range for flux limits
  294 + GEnergy emin = GEnergy(m_emin, "TeV");
  295 + GEnergy emax = GEnergy(m_emax, "TeV");
  296 +
  297 + // Compute upper limit intensity and fluxes
  298 + m_diff_ulimit = m_skymodel->spectral()->eval(eref, GTime());
  299 + m_flux_ulimit = m_skymodel->spectral()->flux(emin, emax);
286 300 m_eflux_ulimit = m_skymodel->spectral()->eflux(emin, emax);
287 301  
288 302 // Write results to logfile
289 303 if (logTerse()) {
290 304 log << std::endl;
291   - log.header1("Upper limit computation finished");
292   - log << "Upper limits are:" << std::endl;
293   - log << " Differential at ";
294   - log << m_eref <<" TeV: ";
295   - log << m_diff_ulimit << " ph/cm2/s/MeV" << std::endl;
296   - log << " Integral [";
297   - log << m_emin << " - " << m_emax;
298   - log << "] TeV: ";
299   - log << m_flux_ulimit << " ph/cm2/s" <<std::endl;
300   - log << " Energy flux [";
301   - log << m_emin << " - " << m_emax;
302   - log << "] TeV: ";
303   - log << m_eflux_ulimit << " erg/cm2/s" <<std::endl;
  305 + log.header1("Upper limit results");
  306 + log << gammalib::parformat("Differential flux limit");
  307 + log << m_diff_ulimit;
  308 + log << " ph/cm2/s/MeV at ";
  309 + log << m_eref << " TeV";
  310 + log << std::endl;
  311 + log << gammalib::parformat("Integral flux limit");
  312 + log << m_flux_ulimit;
  313 + log << " ph/cm2/s within [";
  314 + log << m_emin << "-" << m_emax;
  315 + log << "] TeV";
  316 + log << std::endl;
  317 + log << gammalib::parformat("Energy flux limit");
  318 + log << m_eflux_ulimit;
  319 + log << " erg/cm2/s within [";
  320 + log << m_emin << "-" << m_emax;
  321 + log << "] TeV";
  322 + log << std::endl;
304 323 }
305 324  
306 325 // Recover original models
... ... @@ -310,119 +329,6 @@ void ctulimit::run(void)
310 329 return;
311 330 }
312 331  
313   -/***********************************************************************//**
314   - * @brief Evaluates the likelihood
315   - *
316   - * @param[in] value Factorised parameter value of the parameter of interest
317   - * @return likelihood value
318   - *
319   - * This method evaluates the likelihood as a function of the parameter of interest
320   - ***************************************************************************/
321   -double ctulimit::evaluate(const double& value)
322   -{
323   - // Initialise likelihood value
324   - double LogL = 0.0;
325   -
326   - // Check if given parameter is within boundaries
327   - if (value > (*m_skymodel->spectral())[0].factor_min() && value < (*m_skymodel->spectral())[0].factor_max()) {
328   -
329   - // Remove source model
330   - m_models.remove(m_srcname);
331   -
332   - // Change spectral model value
333   - (*m_skymodel->spectral())[0].factor_value(value);
334   -
335   - // Re-append model to container
336   - m_models.append(*m_skymodel);
337   -
338   - // Assign new models to observations
339   - m_obs.models(m_models);
340   -
341   - // Evaluate likelihood for new model container
342   - m_obs.eval();
343   -
344   - // Retrieve likelihood
345   - LogL = m_obs.logL();
346   -
347   - } // endif: value was inside allowed range
348   -
349   - else {
350   - // throw exception if value out of range
351   - std::string msg = "Value out of range requested: "
352   - "To omit this error, you could increase the allowed "
353   - "parameter range in the model xml file.";
354   - throw GException::invalid_value(G_EVALUATE, msg);
355   - }
356   -
357   - // Compute function value
358   - double difflike = LogL - m_bestloglike - m_dloglike;
359   -
360   - // Return
361   - return difflike;
362   -}
363   -
364   -/***********************************************************************//**
365   - * @brief Induces upper limit computation by using a bisection method
366   - *
367   - * This method calculates the upper limit using a bisection method
368   - ***************************************************************************/
369   -void ctulimit::ulimit_bisection(const double& min, const double& max, const double& scale)
370   -{
371   - // copy values to working values
372   - double wrk_min = min;
373   - double wrk_max = max;
374   -
375   - // Initialise counter
376   - int iter=0;
377   -
378   - // Loop until breaking condition is reached
379   - while (true) {
380   -
381   - // Log information
382   - if (logExplicit()) {
383   - log << " Iteration ";
384   - log << iter;
385   - log << ": Parameter range reduced to [";
386   - log << wrk_min * scale;
387   - log <<"; ";
388   - log << wrk_max * scale;
389   - log << "]" << std::endl;
390   - }
391   -
392   - // Throw exception if maximum iterations are reached
393   - if( iter > m_max_iter) {
394   - throw GException::invalid_value(G_UL_BISECTION, "Maximum iterations reached");;
395   - }
396   -
397   - // compute center of boundary
398   - double mid = (wrk_min + wrk_max) / 2.0;
399   -
400   - // Calculate function value
401   - double eval_mid = evaluate(mid);
402   -
403   - // Check for convergence inside tolerance
404   - if (std::abs(eval_mid) < m_tol) {
405   - break;
406   - }
407   -
408   - // change boundaries for further iteration
409   - if (eval_mid > 0.0) {
410   - wrk_max = mid;
411   - }
412   - else if (eval_mid < 0.0) {
413   - wrk_min = mid;
414   - }
415   -
416   - // increment counter
417   - iter++;
418   -
419   - } // endwhile
420   -
421   - // Return
422   - return;
423   -
424   -}
425   -
426 332  
427 333 /***********************************************************************//**
428 334 * @brief Save maps
... ... @@ -466,26 +372,27 @@ void ctulimit::save(void)
466 372 ***************************************************************************/
467 373 void ctulimit::init_members(void)
468 374 {
469   - // Initialise members
  375 + // Initialise user parameters
470 376 m_outfile.clear();
471 377 m_srcname.clear();
472   - m_bestloglike = 0.0;
473   - m_flux_ulimit = 0.0;
474   - m_diff_ulimit = 0.0;
475   - m_eflux_ulimit = 0.0;
476   - m_dloglike = 0.0;
477   - m_tol = 1e-6;
478   - m_max_iter = 50;
479   - m_eref = 0.0;
480   - m_emin = 0.0;
481   - m_emax = 0.0;
482   - m_sigma_min = 0.0;
483   - m_sigma_max = 0.0;
484   - m_skymodel = NULL;
  378 + m_sigma_min = 0.0;
  379 + m_sigma_max = 0.0;
  380 + m_eref = 0.0;
  381 + m_emin = 0.0;
  382 + m_emax = 0.0;
  383 + m_tol = 1.0e-6;
  384 + m_max_iter = 50;
485 385  
486 386 // Initialise protected members
487 387 m_obs.clear();
488 388 m_models.clear();
  389 + m_dlogL = 0.0;
  390 + m_skymodel = NULL;
  391 + m_model_par = NULL;
  392 + m_best_logL = 0.0;
  393 + m_flux_ulimit = 0.0;
  394 + m_diff_ulimit = 0.0;
  395 + m_eflux_ulimit = 0.0;
489 396  
490 397 // Return
491 398 return;
... ... @@ -499,28 +406,29 @@ void ctulimit::init_members(void)
499 406 ***************************************************************************/
500 407 void ctulimit::copy_members(const ctulimit& app)
501 408 {
502   - // Copy attributes
503   - m_outfile = app.m_outfile;
504   - m_srcname = app.m_srcname;
505   - m_diff_ulimit = app.m_diff_ulimit;
506   - m_flux_ulimit = app.m_flux_ulimit;
507   - m_eflux_ulimit = app.m_eflux_ulimit;
508   - m_bestloglike = app.m_bestloglike;
509   - m_dloglike = app.m_dloglike;
510   - m_tol = app.m_tol;
511   - m_max_iter = app.m_max_iter;
512   - m_eref = app.m_eref;
513   - m_emin = app.m_emin;
514   - m_emax = app.m_emax;
515   - m_sigma_min = app.m_sigma_min;
516   - m_sigma_max = app.m_sigma_max;
  409 + // Copy user parameters
  410 + m_outfile = app.m_outfile;
  411 + m_srcname = app.m_srcname;
  412 + m_sigma_min = app.m_sigma_min;
  413 + m_sigma_max = app.m_sigma_max;
  414 + m_eref = app.m_eref;
  415 + m_emin = app.m_emin;
  416 + m_emax = app.m_emax;
  417 + m_tol = app.m_tol;
  418 + m_max_iter = app.m_max_iter;
  419 +
517 420  
518 421 // Copy protected members
519   - m_obs = app.m_obs;
520   - m_models = app.m_models;
  422 + m_obs = app.m_obs;
  423 + m_models = app.m_models;
  424 + m_dlogL = app.m_dlogL;
  425 + m_best_logL = app.m_best_logL;
  426 + m_diff_ulimit = app.m_diff_ulimit;
  427 + m_flux_ulimit = app.m_flux_ulimit;
  428 + m_eflux_ulimit = app.m_eflux_ulimit;
521 429  
522   - // Clone protected members
523   - m_skymodel = (app.m_skymodel != NULL) ? app.m_skymodel->clone() : NULL;
  430 + // Extract model parameter
  431 + get_model_parameter();
524 432  
525 433 // Return
526 434 return;
... ... @@ -541,14 +449,10 @@ void ctulimit::free_members(void)
541 449 * @brief Get application parameters
542 450 *
543 451 * @exception GException::invalid_value
544   - * Test source not found or no RA/DEC parameters found for test
545   - * source.
  452 + * Test source not found.
546 453 *
547 454 * Get all task parameters from parameter file or (if required) by querying
548   - * the user. Most parameters are only required if no observation exists so
549   - * far in the observation container. In this case, a single CTA observation
550   - * will be added to the container, using the definition provided in the
551   - * parameter file.
  455 + * the user.
552 456 ***************************************************************************/
553 457 void ctulimit::get_parameters(void)
554 458 {
... ... @@ -586,16 +490,20 @@ void ctulimit::get_parameters(void)
586 490 throw GException::invalid_value(G_GET_PARAMETERS, msg);
587 491 }
588 492  
589   - double CL = (*this)["cl"].real();
590   - if (CL != 0.95) {
591   - std::string msg = "Confidence level different from 95% requested"
592   - "Currently only 95% CL is possible";
  493 + // Get relevant model and parameter for upper limit computation
  494 + get_model_parameter();
  495 +
  496 + // Get confidence level
  497 + double confidence = (*this)["confidence"].real();
  498 + if (confidence != 0.95) {
  499 + std::string msg = "Confidence level different from 95% requested."
  500 + " Currently only 95% is supported.";
593 501 throw GException::invalid_value(G_GET_PARAMETERS, msg);
594 502 }
595 503 else {
596 504 // Set Likelihood difference for 95% CL.
597 505 // See Minuit Handbook
598   - m_dloglike = 3.84 / 2.0;
  506 + m_dlogL = 3.84 / 2.0;
599 507 }
600 508  
601 509 // Read starting boundaries for bisection
... ... @@ -608,7 +516,7 @@ void ctulimit::get_parameters(void)
608 516 m_emax = (*this)["emax"].real();
609 517  
610 518 // Read precision
611   - m_tol = (*this)["tol"].real();
  519 + m_tol = (*this)["tol"].real();
612 520 m_max_iter = (*this)["max_iter"].integer();
613 521  
614 522 // Optionally read ahead parameters so that they get correctly
... ... @@ -622,3 +530,156 @@ void ctulimit::get_parameters(void)
622 530 }
623 531  
624 532  
  533 +/***********************************************************************//**
  534 + * @brief Get application parameters
  535 + *
  536 + * @exception GException::invalid_value
  537 + * Did not find a valid model parameter
  538 + *
  539 + * Extracts a pointer to the sky model (m_skymodel) and a pointer to the
  540 + * relevant model parameter (m_model_par) from the model container.
  541 + ***************************************************************************/
  542 +void ctulimit::get_model_parameter(void)
  543 +{
  544 + // Get relevant model and parameter for upper limit computation.
  545 + GModels& models = const_cast<GModels&>(m_obs.models());
  546 + m_skymodel = dynamic_cast<GModelSky*>(models[m_srcname]);
  547 + if (m_skymodel == NULL) {
  548 + std::string msg = "Source \""+m_srcname+"\" is not a sky model. "
  549 + "Please specify the name of a sky model for "
  550 + "upper limit computation.";
  551 + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
  552 + }
  553 + if (m_skymodel->spectral()->has_par("Normalization")) {
  554 + m_model_par = &(m_skymodel->spectral()->operator[]("Normalization"));
  555 + }
  556 + else if (m_skymodel->spectral()->has_par("Prefactor")) {
  557 + m_model_par = &(m_skymodel->spectral()->operator[]("Prefactor"));
  558 + }
  559 + else if (m_skymodel->spectral()->has_par("Integral")) {
  560 + m_model_par = &(m_skymodel->spectral()->operator[]("Integral"));
  561 + }
  562 + else {
  563 + std::string msg = "Require spectral parameter \"Normalization\", "
  564 + "\"Prefactor\" or \"Integral\" for upper limit "
  565 + "computation. The specified source \""+m_srcname+
  566 + "\" does not have such a parameter.";
  567 + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
  568 + }
  569 +
  570 + // Return
  571 + return;
  572 +}
  573 +
  574 +
  575 +/***********************************************************************//**
  576 + * @brief Performs upper limit computation by using a bisection method
  577 + *
  578 + * @param[in] min Minimum parameter value
  579 + * @param[in] max Maximum parameter value
  580 + *
  581 + * This method calculates the upper limit using a bisection method.
  582 + ***************************************************************************/
  583 +void ctulimit::ulimit_bisection(const double& min, const double& max)
  584 +{
  585 + // Copy values to working values
  586 + double wrk_min = min;
  587 + double wrk_max = max;
  588 +
  589 + // Initialise iteration counter
  590 + int iter = 0;
  591 +
  592 + // Loop until breaking condition is reached
  593 + while (true) {
  594 +
  595 + // Log information
  596 + if (logExplicit()) {
  597 + log << gammalib::parformat("Iteration "+gammalib::str(iter));
  598 + log << "[";
  599 + log << wrk_min;
  600 + log << ", ";
  601 + log << wrk_max;
  602 + log << "]" << std::endl;
  603 + }
  604 +
  605 + // Throw exception if maximum iterations are reached
  606 + if (iter > m_max_iter) {
  607 + std::string msg = "The maximum number of "+gammalib::str(m_max_iter)+
  608 + " has been reached. You may consider to increase"
  609 + " the \"max_iter\" parameter and re-run ctulimit.";
  610 + throw GException::invalid_value(G_UL_BISECTION, msg);
  611 + }
  612 +
  613 + // Compute center of boundary
  614 + double mid = (wrk_min + wrk_max) / 2.0;
  615 +
  616 + // Calculate function value
  617 + double eval_mid = evaluate(mid);
  618 +
  619 + // Check for convergence inside tolerance
  620 + if (std::abs(eval_mid) < m_tol) {
  621 + break;
  622 + }
  623 +
  624 + // Change boundaries for further iteration
  625 + if (eval_mid > 0.0) {
  626 + wrk_max = mid;
  627 + }
  628 + else if (eval_mid < 0.0) {
  629 + wrk_min = mid;
  630 + }
  631 +
  632 + // Increment counter
  633 + iter++;
  634 +
  635 + } // endwhile
  636 +
  637 + // Return
  638 + return;
  639 +
  640 +}
  641 +
  642 +
  643 +/***********************************************************************//**
  644 + * @brief Evaluates the log-likelihood
  645 + *
  646 + * @param[in] value Parameter factor value
  647 + * @return Log-likelihood value
  648 + *
  649 + * This method evaluates the log-likelihood as a function of the parameter
  650 + * of interest.
  651 + ***************************************************************************/
  652 +double ctulimit::evaluate(const double& value)
  653 +{
  654 + // Initialise log-likelihood value
  655 + double logL = 0.0;
  656 +
  657 + // Check if given parameter is within boundaries
  658 + if (value > m_model_par->min() && value < m_model_par->max()) {
  659 +
  660 + // Change parameter factor
  661 + m_model_par->value(value);
  662 +
  663 + // Evaluate likelihood for new model container
  664 + m_obs.eval();
  665 +
  666 + // Retrieve likelihood
  667 + logL = m_obs.logL();
  668 +
  669 + } // endif: value was inside allowed range
  670 +
  671 + // ... otherwise signal that the parameter went outside the boundaries
  672 + else {
  673 + std::string msg = "Value of parameter \""+m_model_par->name()+"\" "
  674 + "outside of validity range requested. To omit "
  675 + "this error please enlarge the parameter range "
  676 + "in the model XML file.";
  677 + throw GException::invalid_value(G_EVALUATE, msg);
  678 + }
  679 +
  680 + // Compute function value
  681 + double logL_difference = logL - m_best_logL - m_dlogL;
  682 +
  683 + // Return
  684 + return logL_difference;
  685 +}
... ...
src/ctulimit/ctulimit.hpp
  View file @411b731
1 1 /***************************************************************************
2   - * ctulimit - Upper limit calculation tool *
  2 + * ctulimit - Upper limit calculation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Michael Mayer *
  4 + * copyright (C) 2015 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file ctulimit.hpp
23   - * @brief upper limit calculation tool interface definition
  23 + * @brief Upper limit calculation tool interface definition
24 24 * @author Michael Mayer
25 25 */
26 26  
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTULIMIT_NAME "ctulimit"
37   -#define CTULIMIT_VERSION "00-01-00"
  37 +#define CTULIMIT_VERSION "01-00-00"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -42,14 +42,13 @@
42 42 *
43 43 * @brief Upper limit calculation tool
44 44 *
45   - * This class computes and upper limit for a given parameter.
  45 + * This class computes and upper limit for a given source.
46 46 *
47 47 * The class operates on predefined observation containers, an individual
48 48 * event list or an observation definition XML file.
49 49 *
50   - * During the computation the likelihood function is inspected to find the best
51   - * value for the upper limit.
52   - *
  50 + * During the computation the likelihood function is inspected to find the
  51 + * best value for the upper limit.
53 52 ***************************************************************************/
54 53 class ctulimit : public ctool {
55 54  
... ... @@ -69,43 +68,41 @@ public:
69 68 void run(void);
70 69 void save(void);
71 70 const GObservations& obs(void) const;
72   - const double& diff_ulimit(void) const;
73   - const double& flux_ulimit(void) const;
74   - const double& eflux_ulimit(void) const;
  71 + const double& diff_ulimit(void) const;
  72 + const double& flux_ulimit(void) const;
  73 + const double& eflux_ulimit(void) const;
75 74  
76 75 protected:
77 76 // Protected methods
78   - void init_members(void);
79   - void copy_members(const ctulimit& app);
80   - void free_members(void);
81   - void get_parameters(void);
82   - void ulimit_bisection(const double& parmin, const double& parmax, const double& scale);
  77 + void init_members(void);
  78 + void copy_members(const ctulimit& app);
  79 + void free_members(void);
  80 + void get_parameters(void);
  81 + void get_model_parameter(void);
  82 + void ulimit_bisection(const double& parmin, const double& parmax);
83 83 double evaluate(const double& value);
84 84  
85 85 // User parameters
86   - std::string m_srcname; //!< Name of source which is moved around
87   - std::string m_outfile; //!< Output ascii file
88   -
89   - GModelSky* m_skymodel; //!< Pointer to variable spectral model
90   -
91   - int m_max_iter; //!< maximum number of iterations
92   - double m_dloglike; //!< Likelihood difference for upper limit computation
93   - double m_bestloglike ; //!< Best fit log likelihood of given model
94   - double m_tol; //!< tolerance for limit determination
95   - double m_sigma_min; //!< Starting value minimum (multiple fit errors above fit values)
96   - double m_sigma_max; //!< Starting value maximum (multiple fit errors above fit values)
97   - double m_eref; //!< reference energy for flux limits
98   - double m_emin; //!< minimum energy for flux limits
99   - double m_emax; //!< maximum energy for flux limits
100   -
101   - double m_diff_ulimit; //!<Differential upper limit value
102   - double m_flux_ulimit; //!< Flux upper limit value
103   - double m_eflux_ulimit; //!< Energy flux upper limits
  86 + std::string m_srcname; //!< Name of source which is moved around
  87 + std::string m_outfile; //!< Output ascii file
  88 + double m_sigma_min; //!< Starting value minimum (multiple fit errors above fit values)
  89 + double m_sigma_max; //!< Starting value maximum (multiple fit errors above fit values)
  90 + double m_eref; //!< Reference energy for flux limits (TeV)
  91 + double m_emin; //!< Minimum energy for flux limits (TeV)
  92 + double m_emax; //!< Maximum energy for flux limits (TeV)
  93 + double m_tol; //!< Tolerance for limit determination
  94 + int m_max_iter; //!< Maximum number of iterations
104 95  
105 96 // Protected members
106   - GObservations m_obs; //!< Observation container
107   - GModels m_models;
108   -
  97 + GObservations m_obs; //!< Observation container
  98 + GModels m_models;
  99 + double m_dlogL; //!< Likelihood difference for upper limit computation
  100 + GModelSky* m_skymodel; //!< Pointer to sky model
  101 + GModelPar* m_model_par; //!< Pointer to model parameter
  102 + double m_best_logL; //!< Best fit log likelihood of given model
  103 + double m_diff_ulimit; //!< Differential upper limit value
  104 + double m_flux_ulimit; //!< Flux upper limit value
  105 + double m_eflux_ulimit; //!< Energy flux upper limits
109 106 };
110 107  
111 108  
... ...
src/ctulimit/ctulimit.par
  View file @411b731
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctulimit.par #
8   -# Version: 00-01-00 #
  8 +# Version: 01-00-00 #
9 9 # #
10 10 # Author: Michael Mayer #
11 11 # michael.mayer@physik.hu-berlin.de #
... ... @@ -17,26 +17,26 @@
17 17 #
18 18 # General parameters
19 19 #===================
20   -inobs, f, a, "events.fits",,,"Input event list, counts cube or observation definition file"
  20 +inobs, f, a, "events.fits",,, "Input event list, counts cube or observation definition file"
21 21 inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
22   -srcname, s, a, "Crab",,,"Source of interest"
  22 +srcname, s, a, "Crab",,, "Source of interest"
23 23 expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
24 24 psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
25 25 caldb, s, a, "dummy",,, "Calibration database"
26 26 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
27   -outfile, f, a, "upperlimit.txt",,, "Output ascii file"
  27 +outfile, f, a, "upperlimit.txt",,, "Output ASCII file"
28 28  
29 29 #
30 30 # Upper limit calculation parameters
31 31 #===================================
32   -cl, r, h, 0.95,,,"Confidence Level"
33   -sigma_min,r, h, 0.0,,, "Start minimum value (multiple sigmas above fit value)"
34   -sigma_max,r, h, 10.0,,, "Start maximum value (multiple sigmas above fit value)"
35   -eref, r, h, 1.0,,, "Reference energy for differential limits"
36   -emin, r, h, 1.0,,, "Minimum energy for flux limits"
37   -emax, r, h, 100.0,,, "Maximum energy for flux limits"
38   -tol, r, h, 1e-5,,, "Calculation tolerance"
39   -max_iter, i, h, 50,,, "Maximum number of iterations"
  32 +confidence, r, h, 0.95,0.0,1.0, "Confidence level"
  33 +sigma_min, r, h, 0.0,,, "Start minimum value (multiple sigmas above fit value)"
  34 +sigma_max, r, h, 10.0,,, "Start maximum value (multiple sigmas above fit value)"
  35 +eref, r, h, 1.0,,, "Reference energy for differential limits (TeV)"
  36 +emin, r, h, 1.0,,, "Minimum energy for flux limits (TeV)"
  37 +emax, r, h, 100.0,,, "Maximum energy for flux limits (TeV)"
  38 +tol, r, h, 1e-5,,, "Computation tolerance"
  39 +max_iter, i, h, 50,,, "Maximum number of iterations"
40 40  
41 41 #
42 42 # Standard parameters
... ...
src/ctulimit/main.cpp
  View file @411b731
1 1 /***************************************************************************
2   - * ctulimit - upper limit calculation tool *
  2 + * ctulimit - upper limit calculation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Michael Mayer *
  4 + * copyright (C) 2015 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
test/test_ctulimit.py
  View file @411b731
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This scripts performs unit tests for the cttsmap tool.
  3 +# This scripts performs unit tests for the ctulimit tool.
4 4 #
5   -# Copyright (C) 2014 Juergen Knoedlseder
  5 +# Copyright (C) 2015 Michael Mayer
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ...