Merge branch '1498-ctbutterfly-should-consider-all-correlations' into 1498-fix-ctbutterfly

Conflicts: examples/show_butterfly.py

Merge branch '1498-ctbutterfly-should-consider-all-correlations' into 1498-fix-ctbutterfly
Conflicts: examples/show_butterfly.py
Jürgen Knödlseder
2 parents 42c52ade 1cf7b40b
Showing 3 changed files with 105 additions and 32 deletions
examples/show_butterfly.py
src/ctbutterfly/ctbutterfly.cpp
src/ctbutterfly/ctbutterfly.hpp
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This script shows how to plot a butterfly created with ctbutterfly 
 #
-# Copyright (C) 2014 Michael Mayer
+# Copyright (C) 2014-2015 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -63,14 +63,19 @@ if __name__ == &quot;__main__&quot;:
         index = nrows - 1 - row
         butterfly_x.append(csv.real(index,0))
-        butterfly_y.append(csv.real(index,1)-csv.real(index,2))    
-
+        low_error = csv.real(index,1)-csv.real(index,2)
+        if low_error < 1e-26:
+            low_error = 1e-26
+        butterfly_y.append(low_error)   
+    
     # plot the butterfly and spectral line       
     plt.figure()
+    plt.ylim([1e-26,1e-14])
     plt.loglog()
     plt.grid()
     plt.plot(line_x,line_y,color='black',ls='-')
     plt.fill(butterfly_x,butterfly_y,color='green',alpha=0.5)
+
     plt.xlabel("Energy [MeV]")
     plt.ylabel(r"dN/dE [MeV$^{-1}$ s$^{-1}$ cm$^{-2}$]")    
     plt.show()
@@ -36,6 +36,7 @@
 /* __ Method name definitions ____________________________________________ */
 #define G_GET_PARAMETERS                      "ctbutterfly::get_parameters()"
 #define G_SAVE                                          "ctbutterfly::save()"
+#define G_GET_MODEL_PARAMETER            "ctbutterfly::get_model_parameter()"
 /* __ Debug definitions __________________________________________________ */
@@ -291,47 +292,38 @@ void ctbutterfly::run(void)
             log.header2(msg);
         }
-        // Initialise model flux value
+        // Initialise flux value
         double model_flux = 0.0;
         // Loop over models
         for (int j = 0; j < models.size(); ++j) {
-            // Check wether model is a skymodel
-            GModelSky* skymodel = dynamic_cast<GModelSky*>(models[j]);
+            // If source name is found, retrieve the nominal flux
+            if (models[j]->name() == m_srcname) {
-            // Yes ...
-            if (skymodel != NULL) {
-
-                // Skip spatial models
-                num_gradient += skymodel->spatial()->size();
-
-                // Get pointer to spectral model
+                // Note that we have to compute the gradients here
+                // repetively to update the gradients in the vector
+                // No check for valid pointers necessary here as this 
+                // was covered by get_model_parameter()
+                GModelSky*      skymodel = dynamic_cast<GModelSky*>(models[j]);
                 GModelSpectral* spectral = skymodel->spectral();
+                model_flux = spectral->eval_gradients(energy, time);
-                // Set flux value of source of interest
-                if (skymodel->name() == m_srcname) {
-                    model_flux = spectral->eval_gradients(energy, time);
-                } // endif: model was source of interest
-                
-                // Loop over model parameters, get gradients
-                // and assign them to the vector
-                for (int k = 0; k < spectral->size(); ++k) {
-                    grad[num_gradient] = (*spectral)[k].gradient();
-                    num_gradient++;
-                } // endfor: looped over spectral parameters
+            }
-                // Skip temporal models
-                num_gradient += skymodel->temporal()->size();
+            // Loop over model parameters
+            for (int i = 0; i < models[j]->size(); ++i) {
-            } // endif: model was sky model
+                // Get model parameter
+                GModelPar parameter = models[j]->at(i);
-            else {
+                // Set gradient value to the vector
+                grad[num_gradient] = parameter.gradient();
-                // Skip other models (e.g. GModelData instances)
-                num_gradient += models[j]->size();
+                // increment vector index
+                num_gradient++;
-            } // endelse: Model was not GModelSky
+            } // endfor: Looped over model parameters
         } // endfor: Looped over models
@@ -341,6 +333,9 @@ void ctbutterfly::run(void)
         // Get the error from the scalar product
         double error = std::sqrt(grad * vector);
+        // Multiply in scale of spectral parameter
+        error *= m_par_scale;
+
         // Store flux, value and energy for saving
         m_fluxes.push_back(model_flux);
         m_energies.push_back(energy.MeV());
@@ -425,6 +420,8 @@ void ctbutterfly::init_members(void)
     m_apply_edisp = false;
     // Initialise protected members
+    m_skymodel  = NULL;
+    m_par_scale = 0.0;
     m_obs.clear();
     m_covariance.clear();
     m_energies.clear();
@@ -456,6 +453,9 @@ void ctbutterfly::copy_members(const ctbutterfly&amp; app)
     m_fluxes     = app.m_fluxes;
     m_errors     = app.m_errors;
+    // Extract model parameter
+    get_model_parameter();
+
     // Return
     return;
 }
@@ -517,9 +517,12 @@ void ctbutterfly::get_parameters(void)
         std::string msg = "Source \""+m_srcname+"\" not found in model "
                           "container. Please add a source with that name "
                           "or check for possible typos.";
-    	throw GException::invalid_value(G_GET_PARAMETERS, msg);
+        throw GException::invalid_value(G_GET_PARAMETERS, msg);
     }
+    // Get relevant model and parameter for model computation
+    get_model_parameter();
+
     // Create energy boundaries from user parameters
     m_ebounds = create_ebounds();
@@ -545,3 +548,65 @@ void ctbutterfly::get_parameters(void)
     // Return
     return;
 }
+
+
+/***********************************************************************//**
+ * @brief Get application parameters
+ *
+ * @exception GException::invalid_value
+ *            Did not find a valid model parameter
+ *
+ * Extracts a pointer to the sky model (m_skymodel) and a pointer to the
+ * relevant model parameter (m_model_par) from the model container.
+ ***************************************************************************/
+void ctbutterfly::get_model_parameter(void)
+{
+    // Continue only if source model exists
+    if (m_obs.models().contains(m_srcname)) {
+
+        // Get relevant model and parameter for upper limit computation.
+        GModels& models = const_cast<GModels&>(m_obs.models());
+        m_skymodel      = dynamic_cast<GModelSky*>(models[m_srcname]);
+        if (m_skymodel == NULL) {
+            std::string msg = "Source \""+m_srcname+"\" is not a sky model. "
+                              "Please specify the name of a sky model for "
+                              "butterfly computation.";
+            throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
+        }
+
+        GModelPar* modelpar;
+
+        if (m_skymodel->spectral()->type() == "NodeFunction") {
+            std::string msg = "\"NodeFunction\" cannot be used as spectral "
+                              "model for an butterfly computation. "
+                              "Please specify another spectral model.";
+            throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
+        }
+        if (m_skymodel->spectral()->has_par("Normalization")) {
+            modelpar = &(m_skymodel->spectral()->operator[]("Normalization"));
+        }
+        else if (m_skymodel->spectral()->has_par("Prefactor")) {
+            modelpar = &(m_skymodel->spectral()->operator[]("Prefactor"));
+        }
+        else if (m_skymodel->spectral()->has_par("Integral")) {
+            modelpar = &(m_skymodel->spectral()->operator[]("Integral"));
+        }
+        else if (m_skymodel->spectral()->has_par("Value")) {
+            modelpar = &(m_skymodel->spectral()->operator[]("Value"));
+        }
+        else {
+            std::string msg = "Require spectral parameter \"Normalization\", "
+                              "\"Prefactor\", \"Integral\" or \"Value\" for "
+                              "butterfly computation. The specified source "
+                              "\""+m_srcname+"\" does not have such a parameter.";
+            throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
+        }
+
+        // Set scale value from parameter
+        m_par_scale = modelpar->scale();
+
+    } // endif: source model existed
+
+    // Return
+    return;
+}
@@ -78,6 +78,7 @@ protected:
     void copy_members(const ctbutterfly& app);
     void free_members(void);
     void get_parameters(void);
+    void get_model_parameter(void);
     // User parameters
     std::string m_srcname;     //!< Name of source to compute butterfly
@@ -88,6 +89,8 @@ protected:
     // Protected members
     GObservations       m_obs;        //!< Observation container
     GMatrixSparse       m_covariance; //!< Covariance matrix
+    GModelSky*          m_skymodel;   //!< Pointer to sky model
+    double              m_par_scale;  //!< Scale of normalising parameter
     std::vector<double> m_energies;   //!< Energy values for storage
     std::vector<double> m_fluxes;     //!< Flux values per energy bin
     std::vector<double> m_errors;     //!< Flux errors per energy bin