diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp index 6db4459..d01661e 100644 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -26,54 +26,36 @@ /***********************************************************************//** - * @brief Setup Crab point source powerlaw model. - ***************************************************************************/ -GModels crab_plaw(void) -{ - // Setup GModels for optimizing - GModelSpatialPtsrc point_source; - GModelSpectralPlaw power_law; - GModel crab; - GModels models; - try { - GSkyDir dir; - //dir.radec_deg(83.6331, +22.0145); - dir.radec_deg(117.0, -33.0); // Adapt to source position in file - point_source = GModelSpatialPtsrc(dir); - power_law = GModelSpectralPlaw(1.0e-7, -2.1); - power_law.par(0)->min(1.0e-12); - crab = GModel(point_source, power_law); - crab.name("Crab"); - models.append(crab); - } - catch (std::exception &e) { - std::cout << std::endl - << "TEST ERROR: Unable setup GModels for optimizing." - << std::endl; - std::cout << e.what() << std::endl; - throw; - } - - // Return model - return models; -} - - -/***********************************************************************//** * @brief Run maximum likelihood analysis ***************************************************************************/ int ctlike::run(void) { // Initialise return code int rc = 0; - + + // Get standard parameters + m_method = toupper(par("method")->value()); + m_caldb = par("caldb")->value(); + m_irf = par("irf")->value(); + m_srcmdl = par("srcmdl")->value(); + m_models = GModels(m_srcmdl); + + // Set fixed parameters + m_max_iter = 1000; + m_opt = new GOptimizerLM; + + // Set optimizer parameters + ((GOptimizerLM*)m_opt)->max_iter(m_max_iter); + // Branch on method - std::string method = par("method")->value(); - if (toupper(method) == "UNBINNED") + if (m_method == "UNBINNED") rc = unbinned(); else rc = binned(); - + + // Free optimizer + delete m_opt; + // Return return rc; } @@ -89,8 +71,6 @@ int ctlike::unbinned(void) // Get parameters std::string evfile = par("evfile")->value(); - std::string caldb = par("caldb")->value(); - std::string irf = par("irf")->value(); // Declare observations GObservations obs; @@ -118,20 +98,17 @@ int ctlike::unbinned(void) // Load data and response and set ROI, energy range and time range // for analysis run.load_unbinned(evfile); - run.response(irf,caldb); + run.response(m_irf, m_caldb); run.roi(roi); run.gti()->add(tstart, tstop); run.ebounds()->append(emin, emax); obs.append(run); - // Setup GModels for optimizing - GModels models = crab_plaw(); - obs.models(models); + // Setup models for optimizing + obs.models(m_models); // Perform LM optimization - GOptimizerLM opt; - opt.max_iter(1000); - obs.optimize(opt); + obs.optimize(*m_opt); std::cout << obs << std::endl; // Return @@ -149,8 +126,6 @@ int ctlike::binned(void) // Get parameters std::string cntmap = par("cntmap")->value(); - std::string caldb = par("caldb")->value(); - std::string irf = par("irf")->value(); // Declare observations GObservations obs; @@ -158,17 +133,14 @@ int ctlike::binned(void) // Load binned CTA observation run.load_binned(cntmap); - run.response(irf,caldb); + run.response(m_irf, m_caldb); obs.append(run); - // Setup GModels for optimizing - GModels models = crab_plaw(); - obs.models(models); + // Setup models for optimizing + obs.models(m_models); // Perform LM optimization - GOptimizerLM opt; - opt.max_iter(1000); - obs.optimize(opt); + obs.optimize(*m_opt); std::cout << obs << std::endl; // Return diff --git a/src/ctlike/ctlike.hpp b/src/ctlike/ctlike.hpp index 19987e0..c016d58 100644 --- a/src/ctlike/ctlike.hpp +++ b/src/ctlike/ctlike.hpp @@ -20,6 +20,7 @@ #define CTLIKE_HPP /* __ Includes ___________________________________________________________ */ +#include #include "GammaLib.hpp" #include "GCTALib.hpp" @@ -45,6 +46,18 @@ public: int run(void); int unbinned(void); int binned(void); + +protected: + // User parameters + std::string m_method; //!< Likelihood fitting method (binned/unbinned) + std::string m_caldb; //!< Calibration database + std::string m_irf; //!< Instrument response functions + std::string m_srcmdl; //!< Source model XML file + GModels m_models; //!< Source models + + // Fixed parameters + int m_max_iter; //!< Maximum number of iterations + GOptimizer* m_opt; //!< Optimizer }; #endif /* CTLIKE_HPP */ diff --git a/src/ctlike/ctlike.par b/src/ctlike/ctlike.par index cad202e..a6aba17 100644 --- a/src/ctlike/ctlike.par +++ b/src/ctlike/ctlike.par @@ -19,17 +19,18 @@ # # Unbinned analysis #================== -evfile, f, a, "../gammalib/inst/cta/test/data/run_00006028_eventlist_reco.fits.gz",,,"Event data file name" +evfile, f, a, "test/data/run_00006028_eventlist_reco.fits.gz",,,"Event data file name" # # Binned analysis #================== -cntmap, f, a, "../gammalib/inst/cta/test/data/run_00006028_cntmap.fits.gz",,, "Counts map" +cntmap, f, a, "test/data/run_00006028_cntmap.fits.gz",,, "Counts map" # # General parameters #=================== method, s, a, "UNBINNED",,, "Maximum likelihood method" -caldb, s, a, "../gammalib/inst/cta/test/irf",,, "Calibration database" +caldb, s, a, "test/irf",,, "Calibration database" irf, s, a, "kb_E_50h_v3",,, "Instrument response function" +srcmdl, f, a, "test/data/source1.xml",,, "Source model" # # Standard parameters #====================