34 #include "GOptimizer.hpp"
37 #define G_GET_MODEL_PARAMETER "ctulimit::get_model_parameter()"
38 #define G_GET_PARAMETER_BRACKETS "ctulimit::get_parameter_brackets(double&, "\
40 #define G_UL_BISECTION "ctulimit::ul_bisection(double&, double&)"
201 this->GApplication::clear();
235 GModels models_orig =
m_obs.models();
238 GOptimizerLM best_opt;
246 log_header1(TERSE,
"Compute best-fit likelihood");
261 log_value(NORMAL,
"Maximum log likelihood",gammalib::str(
m_best_logL,3));
270 log_header1(TERSE,
"Extract best-fit likelihood");
273 log_value(NORMAL,
"Maximum log likelihood",gammalib::str(
m_best_logL,3));
284 log_header1(TERSE,
"Compute upper limit");
287 log_value(NORMAL,
"Model name",
m_skymodel->name());
288 log_value(NORMAL,
"Parameter name",
m_model_par->name());
290 log_value(NORMAL,
"Confidence level", gammalib::str(
m_confidence*100.0)+
"%");
291 log_value(NORMAL,
"Maximum log likelihood",gammalib::str(
m_best_logL,3));
292 log_value(NORMAL,
"Log-likelihood difference",
m_dlogL);
293 log_value(NORMAL,
"Initial factor value range",
294 "["+gammalib::str(parmin)+
", "+gammalib::str(parmax)+
"]");
300 log_value(NORMAL,
"Upper limit factor value",
m_model_par->factor_value());
301 log_value(NORMAL,
"Upper limit param. value",
m_model_par->value());
307 log_header1(TERSE,
"Upper limit results");
310 log_value(TERSE,
"Differential flux limit",
312 gammalib::str(
m_eref)+
" TeV");
313 log_value(TERSE,
"Integral flux limit",
315 gammalib::str(
m_emin)+
"-"+
316 gammalib::str(
m_emax)+
"] TeV");
317 log_value(TERSE,
"Energy flux limit",
319 gammalib::str(
m_emin)+
"-"+
320 gammalib::str(
m_emax)+
"] TeV");
323 m_obs.models(models_orig);
469 double sigma = gammalib::erfinv(
m_confidence) * gammalib::sqrt_two;
477 m_eref = (*this)[
"eref"].real();
478 m_emin = (*this)[
"emin"].real();
479 m_emax = (*this)[
"emax"].real();
482 m_opt.eps((*
this)[
"like_accuracy"].real());
483 m_opt.max_iter((*
this)[
"max_iter"].integer());
486 m_tol = (*this)[
"tol"].real();
490 m_chatter =
static_cast<GChatter
>((*this)[
"chatter"].integer());
514 const char* pars[] = {
"Normalization",
"Value",
"Prefactor",
515 "Integral",
"PhotonFlux",
"EnergyFlux",
526 std::string msg =
"Source \""+
m_srcname+
"\" not found in model "
527 "container. Please add a source with that name "
528 "or check for possible typos.";
535 GModels& models =
const_cast<GModels&
>(
m_obs.models());
538 std::string msg =
"Source \""+
m_srcname+
"\" is not a sky model. Please "
539 "specify the name of a sky model for upper limit "
554 std::string msg =
"Spatial or spectral model of source \""+
556 "\". Please specify a valid parameter name "
557 "or leave the parameter name blank for "
558 "autodetermination of an intensity- or "
559 "flux-like parameter.";
569 if (
m_skymodel->spectral()->type() ==
"NodeFunction") {
570 std::string msg =
"Spectral model of source \""+
m_srcname+
"\" is "
571 "a \"NodeFunction\" but no parameter name was "
572 "specified. Please use the \"parname\" parameter "
573 "to specify the parameter that should be used "
574 "for the upper limit computation.";
580 for (
const char** par = pars; *par != NULL; ++par) {
591 std::string msg =
"Require one of the following spectral "
592 "parameters for upper limit computation: ";
593 for (
const char** par = pars; *par != NULL; ++par) {
597 msg.append(
"\""+std::string(*par)+
"\"");
599 msg.append(
". The specified source \""+
m_srcname+
"\" does not "
600 "have such a parameter.");
623 log_header1(TERSE,
"Find parameter brackets");
652 std::string msg =
"The maximum number of "+gammalib::str(
m_max_iter)+
653 " has been reached. You may consider to increase"
654 " the \"max_iter\" parameter and re-run ctulimit.";
669 log_value(NORMAL,
"Parameter range",
670 "["+gammalib::str(parmin)+
", "+gammalib::str(parmax)+
"] "
671 "logL("+gammalib::str(parmax)+
")="+gammalib::str(logL));
675 if (eval_mid < 0.0) {
712 double wrk_min = parmin;
713 double wrk_max = parmax;
717 bool restart =
false;
724 std::string msg =
"The maximum number of "+gammalib::str(
m_max_iter)+
725 " has been reached. You may consider to increase"
726 " the \"max_iter\" parameter and re-run ctulimit.";
737 double mid = (wrk_min + wrk_max) / 2.0;
744 log_value(NORMAL,
"Iteration "+gammalib::str(iter),
745 "["+gammalib::str(wrk_min)+
", "+gammalib::str(wrk_max)+
"] "
746 "logL("+gammalib::str(mid)+
")="+gammalib::str(logL)+
" "
747 "dlogL="+gammalib::str(eval_mid));
750 if (std::abs(eval_mid) <
m_tol) {
759 if (std::abs(eps) < 1.0e-6) {
761 " *** WARNING: Parameter range mid-point "+gammalib::str(mid)+
" is "
762 "close to best factor\n"
764 "log-likelihood value "+gammalib::str(logL)+
"\n"
765 " at the mid-point differs significantly from the "
766 "log-likelihood value\n"
767 " "+gammalib::str(
m_best_logL)+
" for the best factor "
769 log_string(TERSE, msg);
774 eps = (mid != 0.0) ? (wrk_min - wrk_max)/mid : wrk_min - wrk_max;
775 if (std::abs(eps) < 1.0e-6) {
777 " *** WARNING: Parameter range ["+gammalib::str(wrk_min)+
", "+
778 gammalib::str(wrk_max)+
"] has reduced\n"
779 " to a narrow interval without reaching convergence! "
781 " log-likelihood value is probably not an absolute "
784 log_string(TERSE, msg);
789 if (!restart && status != 0) {
800 " Adopt "+gammalib::str(logL)+
" as best log-likelihood and restart "
802 log_string(TERSE, msg);
810 if (eval_mid > 0.0) {
842 GEnergy eref = GEnergy(
m_eref,
"TeV");
845 GEnergy emin = GEnergy(
m_emin,
"TeV");
846 GEnergy emax = GEnergy(
m_emax,
"TeV");
850 GModelSpectral* spectral =
m_skymodel->spectral();
851 GModelSpectralNodes* nodes = NULL;
856 GModelSpatialDiffuseCube* cube =
857 dynamic_cast<GModelSpatialDiffuseCube*
>(
m_skymodel->spatial());
862 cube->mc_cone(GSkyRegionCircle(0.0, 0.0, 180.0));
865 nodes =
new GModelSpectralNodes(cube->spectrum());
866 for (
int i = 0; i < nodes->nodes(); ++i) {
867 GEnergy energy = nodes->energy(i);
868 double intensity = nodes->intensity(i);
869 double value = spectral->eval(energy);
870 nodes->intensity(i, value*intensity);
876 if (nodes->nodes() > 0) {
886 if (spectral != NULL) {
#define G_GET_PARAMETER_BRACKETS
Upper limit calculation tool interface implementation.
ctulimit(void)
Void constructor.
GChatter m_chatter
Chattiness.
ctlikelihood & operator=(const ctlikelihood &app)
Assignment operator.
void get_model_parameter(void)
Get model parameter.
double evaluate(GModelPar &par, const double &value)
Evaluates the log-likelihood function.
double m_dlogL
Likelihood difference for upper limit computation.
double m_diff_ulimit
Differential upper limit value.
Upper limit calculation tool.
void copy_members(const ctulimit &app)
Copy class members.
virtual ~ctulimit(void)
Destructor.
double m_emax
Maximum energy for flux limits (TeV)
void set_obs_statistic(const std::string &statistic)
Set fit statistic for CTA observations.
void free_members(void)
Delete class members.
double m_sigma_max
Starting value maximum (multiple fit errors above fit values)
double m_sigma_min
Starting value minimum (multiple fit errors above fit values)
GModelSky * m_skymodel
Pointer to sky model.
std::string m_parname
Name of parameter for upper limit computation.
ctulimit & operator=(const ctulimit &app)
Assignment operator.
bool m_is_spatial
Signal that model parameter is spatial parameter.
void clear(void)
Clear ctulimit tool.
#define G_GET_MODEL_PARAMETER
void free_members(void)
Delete class members.
void free_members(void)
Delete class members.
double m_eflux_ulimit
Energy flux upper limits.
GModelPar * m_model_par
Pointer to model parameter.
void get_parameter_brackets(double &parmin, double &parmax)
Determine parameter brackets.
double m_emin
Minimum energy for flux limits (TeV)
bool m_apply_edisp
Apply energy dispersion?
void init_members(void)
Initialise class members.
void get_parameters(void)
Get application parameters.
double m_confidence
Confidence level.
int m_max_iter
Maximum number of iterations.
void save(void)
Save upper limits.
void init_members(void)
Initialise class members.
double m_best_value
Best parameter value factor.
void init_members(void)
Initialise class members.
double m_tol
Tolerance for limit determination.
void process(void)
Compute upper limit.
GObservations m_obs
Observation container.
double m_eref
Reference energy for flux limits (TeV)
double m_best_error
Best parameter value error.
Base class for likelihood tools.
void compute_ulimit(void)
Compute upper limit intensity and fluxes.
std::string m_srcname
Name of source for upper limit computation.
double m_flux_ulimit
Flux upper limit value.
GOptimizerLM m_opt
Optimizer.
double m_best_logL
Best fit log likelihood of given model.
void ulimit_bisection(const double &parmin, const double &parmax)
Calculate upper limit using a bisection method.