ctbutterfly now takes into account all correlations with other models

Michael Mayer
1 parent 96602a11
Showing 3 changed files with 103 additions and 30 deletions
examples/show_butterfly.py
src/ctbutterfly/ctbutterfly.cpp
src/ctbutterfly/ctbutterfly.hpp
@@ -63,14 +63,20 @@ 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()
 \ No newline at end of file
@@ -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 __________________________________________________ */
@@ -279,47 +280,37 @@ 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 (models[j]->name() == m_srcname) {
-            // Yes ...
-            if (skymodel != NULL) {
-
-                // Skip spatial models
-                num_gradient += skymodel->spatial()->size();
-
-                // Get pointer to spectral model
+                // If source name is found retrieve the nominal flux
+                // 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
@@ -329,6 +320,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());
@@ -404,6 +398,8 @@ void ctbutterfly::init_members(void)
     m_ebounds.clear();
     // Initialise protected members
+    m_skymodel     = NULL;
+    m_par_scale    = 0.0;
     m_obs.clear();
     m_covariance.clear();
     m_energies.clear();
@@ -434,6 +430,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;
 }
@@ -495,9 +494,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();
@@ -520,3 +522,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
@@ -87,6 +88,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