Merge branch 'integration' into devel

Jürgen Knödlseder
2 parents 1b2f2d14 16092b5c
Showing 3 changed files with 46 additions and 29 deletions
pyext/ctools/obsutils.py
scripts/cspull.py
src/cterror/cterror.cpp
@@ -57,7 +57,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
         sim["debug"].boolean(True)
     # Set chatter level
-    sim["chatter"].chatter = chatter
+    sim["chatter"].integer(chatter)
     # Run ctobssim application. This will loop over all observations in the
     # container and simulation the events for each observation. Note that
@@ -105,7 +105,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
             bin["debug"].boolean(True)
         # Set chatter level
-        bin["chatter"].chatter = chatter
+        bin["chatter"].integer(chatter)
         # Run ctbin application. This will loop over all observations in
         # the container and bin the events in counts maps
@@ -157,7 +157,7 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False):
         like["debug"].boolean(True)
     # Set chatter level
-    like["chatter"].chatter = chatter
+    like["chatter"].integer(chatter)
     # Optionally apply energy dispersion
     like["edisp"].boolean(edisp)
@@ -172,12 +172,13 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False):
 # ============================================================== #
 # Fit observations and determine errors using likelihood profile #
 # ============================================================== #
-def cterror(obs, log=False, debug=False, chatter=2):
+def cterror(obs, srcname, log=False, debug=False, chatter=2):
     """
     Perform maximum likelihood fitting of observations in the container.
     Parameters:
-     obs   - Observation container
+     obs     - Observation container
+     srcname - Source name
     Keywords:
      log     - Create log file(s)
      debug   - Create screen dump
@@ -186,6 +187,9 @@ def cterror(obs, log=False, debug=False, chatter=2):
     """
     # Allocate cterror application
     error = ctools.cterror(obs)
+
+    # Set cterror parameters
+    error["srcname"].string(srcname)
     # Optionally open the log file
     if log:
@@ -196,7 +200,7 @@ def cterror(obs, log=False, debug=False, chatter=2):
         error["debug"].boolean(True)
     # Set chatter level
-    error["chatter"].chatter = chatter
+    error["chatter"].integer(chatter)
     # Run cterror application.
     error.run()
@@ -166,8 +166,9 @@ class cspull(ctools.cscript):
         self.m_profile = self["profile"].boolean()
         # Set some fixed parameters
-        self.m_log   = False                    # Logging in client tools
-        self.m_debug = self["debug"].boolean()  # Debugging in client tools
+        self.m_log     = False                     # Logging in client tools
+        self.m_chatter = self["chatter"].integer() # Chatter level 
+        self.m_debug   = self["debug"].boolean()   # Debugging in client tools
         # Return
         return
@@ -276,7 +277,7 @@ class cspull(ctools.cscript):
          seed - Random number generator seed
         """
         # Write header
-        if self.logExplicit():
+        if self.logNormal():
             self.log.header2("Trial "+str(seed+1))
         # Simulate events
@@ -286,7 +287,9 @@ class cspull(ctools.cscript):
                            binsz=self.m_binsz, \
                            npix=self.m_npix, \
                            edisp=self.m_edisp, \
-                           log=self.m_log, debug=self.m_debug)
+                           log=self.m_log, \
+                           debug=self.m_debug, \
+                           chatter=self.m_chatter)
         # Determine number of events in simulation
         nevents = 0.0
@@ -294,7 +297,7 @@ class cspull(ctools.cscript):
             nevents += run.events().number()
         # Write simulation results
-        if self.logExplicit():
+        if self.logNormal():
             self.log.header3("Simulation")
             self.log.parformat("Number of simulated events")
             self.log(nevents)
@@ -302,11 +305,18 @@ class cspull(ctools.cscript):
         # Fit model
         if self.m_profile:
-            like  = obsutils.cterror(obs, \
-                                     log=self.m_log, debug=self.m_debug)
+            models = obs.models()
+            for i in range(models.size()):
+                model_name = models[i].name()
+                like       = obsutils.cterror(obs, model_name, \
+                                              log=self.m_log, \
+                                              debug=self.m_debug,
+                                              chatter=self.m_chatter)
         else:
             like = obsutils.fit(obs, edisp=self.m_edisp, \
-                                log=self.m_log, debug=self.m_debug)
+                                log=self.m_log, \
+                                debug=self.m_debug, \
+                                chatter=self.m_chatter)
         # Store results
         logL   = like.opt().value()
@@ -314,7 +324,7 @@ class cspull(ctools.cscript):
         models = like.obs().models()
         # Write result header
-        if self.logExplicit():
+        if self.logNormal():
             self.log.header3("Pulls")
         # Gather results
@@ -356,7 +366,7 @@ class cspull(ctools.cscript):
                     values["Unc_"+name] = error
                     # Write result
-                    if self.logExplicit():
+                    if self.logNormal():
                         self.log.parformat(name)
                         self.log(pull)
                         self.log(" (")
@@ -37,7 +37,6 @@
 #define G_GET_PARAMETERS                         "cterror::get_parameters()"
 #define G_ERR_BISECTION         "cterror::error_bisection(double&, double&)"
 #define G_EVALUATE                              "cterror::evaluate(double&)"
-#define G_RUN                                               "cterror::run()"
 /* __ Debug definitions __________________________________________________ */
@@ -241,6 +240,9 @@ void cterror::run(void)
     m_obs.errors(*m_opt);
     m_best_logL = m_obs.logL();
+    // Store optimizer for later recovery
+    GOptimizerLM best_opt = *m_opt;
+
     // Write optimised model into logger
     if (logTerse()) {
         log << *m_opt << std::endl;
@@ -255,14 +257,8 @@ void cterror::run(void)
         // Save best fitting models
         GModels models_best = m_obs.models();
-        // Get pointer on sky model
-        GModelSky* model  = dynamic_cast<GModelSky*>(models_best[m_srcname]);
-        if (model == NULL) {
-            std::string msg = "Source \""+m_srcname+"\" is not a sky model. "
-                              "Please specify the name of a sky model for "
-                              "parameter error computation.";
-            throw GException::invalid_value(G_RUN, msg);
-        }
+        // Get pointer on model
+        GModel* model = models_best[m_srcname];
         // Get number of parameters
         int npars = model->size();
@@ -278,10 +274,9 @@ void cterror::run(void)
             // Initialise with best fitting models
             m_obs.models(models_best);
-            // Get pointer on sky model parameter
-            GModels&   current_models = const_cast<GModels&>(m_obs.models());
-            GModelSky* current_model  = dynamic_cast<GModelSky*>(current_models[m_srcname]);
-            m_model_par               = &(current_model->at(i));
+            // Get pointer on model parameter
+            GModels& current_models = const_cast<GModels&>(m_obs.models());
+            m_model_par             = &(current_models[m_srcname]->at(i));
             // Extract current value
             m_value = m_model_par->factor_value();
@@ -360,6 +355,9 @@ void cterror::run(void)
     } // endif: source model exists
+    // Recover optimizer
+    *m_opt = best_opt;
+
     // Return
     return;
 }
@@ -595,6 +593,11 @@ double cterror::error_bisection(const double&amp; min, const double&amp; max)
             break;
         }
+        // Check if interval is smaller than 1.0e-6
+        if (std::abs(wrk_max-wrk_min) < 1.0e-6) {
+            break;
+        }
+
         // If we are on the crescent side of the parabola ...
         if (mid > m_value) {