diff --git a/examples/show_butterfly.py b/examples/show_butterfly.py index ad2774b..ffdabf6 100755 --- a/examples/show_butterfly.py +++ b/examples/show_butterfly.py @@ -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__ == "__main__": 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() diff --git a/src/ctbutterfly/ctbutterfly.cpp b/src/ctbutterfly/ctbutterfly.cpp index 21aaf21..588563a 100644 --- a/src/ctbutterfly/ctbutterfly.cpp +++ b/src/ctbutterfly/ctbutterfly.cpp @@ -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(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(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& 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(m_obs.models()); + m_skymodel = dynamic_cast(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; +} diff --git a/src/ctbutterfly/ctbutterfly.hpp b/src/ctbutterfly/ctbutterfly.hpp index d30617f..c23ba12 100644 --- a/src/ctbutterfly/ctbutterfly.hpp +++ b/src/ctbutterfly/ctbutterfly.hpp @@ -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 m_energies; //!< Energy values for storage std::vector m_fluxes; //!< Flux values per energy bin std::vector m_errors; //!< Flux errors per energy bin