diff --git a/doc/source/reference_manual/ctbutterfly.rst b/doc/source/reference_manual/ctbutterfly.rst
index 5e2e591..0a6b7ca 100644
--- a/doc/source/reference_manual/ctbutterfly.rst
+++ b/doc/source/reference_manual/ctbutterfly.rst
@@ -9,8 +9,15 @@ Computes butterfly.
Synopsis
--------
-Calculates the confidence band of a spectral model, taking into account the
-covariance matrix of a likelihood fit.
+Calculates the confidence band of a power law model, taking into account the
+covariance matrix of a likelihood fit. The tool derives the envelope of
+power law models with parameters that lie within the error ellipse. It
+outputs an ASCII file in column separated value (CSV) format containing
+the results. Each row in the result file corresponds to an energy. The
+columns give the energy in MeV, the fitted intensity in ph/cm2/s/MeV,
+and the minimum and maximum intensity that envelope the power law model
+with a given confidence (by default a confidence level of 68% is used, but
+the level can be adjusted using the ``confidence`` parameter).
General parameters
@@ -45,7 +52,13 @@ General parameters
``outfile [file]``
Output butterfly ASCII file name.
+
+``(refit = no) [boolean]``
+ Performs refitting of solution ignoring any provided covariance matrix.
+``(confidence = 0.68) [real]``
+ Confidence level for error computation.
+
``(matrix = "NONE") [file]``
Input covariance matrix file (not used)
@@ -89,11 +102,12 @@ Standard parameters
``(mode = ql) [string]``
Mode of automatic parameters (default is "ql", i.e. "query and learn").
-``(logfile = ctskymap.log) [string]``
+``(logfile = ctbutterfly.log) [string]``
Name of log file.
Related tools
-------------
-None
+:ref:`ctulimit`
+:ref:`cterror`
diff --git a/doc/source/reference_manual/cterror.rst b/doc/source/reference_manual/cterror.rst
index e57f019..91ec9f1 100644
--- a/doc/source/reference_manual/cterror.rst
+++ b/doc/source/reference_manual/cterror.rst
@@ -89,4 +89,5 @@ Standard parameters
Related tools
-------------
-:ref:`ctulimit`
\ No newline at end of file
+:ref:`ctulimit`
+:ref:`ctbutterfly`
diff --git a/doc/source/reference_manual/ctlike.rst b/doc/source/reference_manual/ctlike.rst
index 348cad2..ccbec47 100644
--- a/doc/source/reference_manual/ctlike.rst
+++ b/doc/source/reference_manual/ctlike.rst
@@ -12,11 +12,10 @@ Synopsis
Determines source model parameters, such as flux, spectral index, source
position, and source extent from a maximum likelihood analysis of CTA data.
The analysis can be done in a binned or an unbinned formulation of the
-log-likelihood function.
-For binned analysis, a counts cube produced by :doc:`ctbin` is required.
-For unbinned analysis, an event list processed by :doc:`ctselect` is required.
-Based on the input file format, ctlike automatically selects between binned
-and unbinned maximum likelihood analysis.
+log-likelihood function. For binned analysis, a counts cube produced by
+:doc:`ctbin` is required. For unbinned analysis, an event list processed by
+:doc:`ctselect` is required. Based on the input file format, ctlike
+automatically selects between binned and unbinned maximum likelihood analysis.
General parameters
diff --git a/examples/show_butterfly.py b/examples/show_butterfly.py
index ffdabf6..ceffc3e 100755
--- a/examples/show_butterfly.py
+++ b/examples/show_butterfly.py
@@ -18,7 +18,6 @@
# along with this program. If not, see .
#
# ==========================================================================
-
import gammalib
import sys
try:
@@ -26,32 +25,34 @@ try:
except:
sys.exit("This script needs matplotlib")
-# =============== #
-# Script entry #
-# =============== #
-if __name__ == "__main__":
- # check for given butterfly file
+# ============ #
+# Script entry #
+# ============ #
+if __name__ == "__main__":
+ """
+ """
+ # Check for given butterfly file
if not len(sys.argv) == 2:
sys.exit("Usage: show_butterfly.py butterfly.txt")
# Read given butterfly file
filename = sys.argv[1]
- csv = gammalib.GCsv(filename)
+ csv = gammalib.GCsv(filename)
- # initialise arrays to be filled
+ # Initialise arrays to be filled
butterfly_x = []
butterfly_y = []
+ line_x = []
+ line_y = []
- line_x = []
- line_y = []
-
- # loop over rows of the file
+ # Loop over rows of the file
nrows = csv.nrows()
for row in range(nrows):
+
# Compute upper edge of confidence band
butterfly_x.append(csv.real(row,0))
- butterfly_y.append(csv.real(row,1)+csv.real(row,2))
+ butterfly_y.append(csv.real(row,2))
# Set line values
line_x.append(csv.real(row,0))
@@ -61,21 +62,19 @@ if __name__ == "__main__":
# of the confidence band
for row in range(nrows):
index = nrows - 1 - row
-
butterfly_x.append(csv.real(index,0))
- low_error = csv.real(index,1)-csv.real(index,2)
+ low_error = csv.real(index,3)
if low_error < 1e-26:
low_error = 1e-26
butterfly_y.append(low_error)
- # plot the butterfly and spectral line
+ # 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.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 588563a..8085375 100644
--- a/src/ctbutterfly/ctbutterfly.cpp
+++ b/src/ctbutterfly/ctbutterfly.cpp
@@ -36,7 +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()"
+#define G_CHECK_MODEL "ctbutterfly::check_model()"
/* __ Debug definitions __________________________________________________ */
@@ -241,11 +241,40 @@ void ctbutterfly::run(void)
log << m_obs << std::endl;
}
+ // If refit is selected then do maximum likelihood fit
+ if (m_refit) {
+
+ // Write header
+ if (logTerse()) {
+ log << std::endl;
+ log.header1("Compute best-fit likelihood");
+ }
+
+ // Optimize and compute errors
+ m_obs.optimize(m_opt);
+ m_obs.errors(m_opt);
+
+ // Get covariance matrix by inverting the curvature matrix. Note that
+ // the curvature matrix is for the unscaled parameters, hence we have
+ // to scale them before we can proceed.
+ GObservations::likelihood likelihood = m_obs.function();
+ m_covariance = likelihood.curvature()->invert();
+
+ // Write optimised model and covariance matrix into logger
+ if (logTerse()) {
+ log << m_opt << std::endl;
+ log << m_obs.models() << std::endl;
+ log << m_covariance << std::endl;
+ }
+
+ } // endif: refit selected
+
// Get model instance for further computations
GModels models = m_obs.models();
- // Compute covariance matrix if necessary
- if (m_covariance.size() == 0) {
+ // If no refit is selected and no covariance matrix has been specified
+ // then compute now the covariance matrix
+ if (!m_refit && m_covariance.size() == 0) {
// Write header
if (logTerse()) {
@@ -258,101 +287,195 @@ void ctbutterfly::run(void)
GObservations::likelihood likelihood = m_obs.function();
likelihood.eval(pars);
- // Get covariance matrix by inverting the curvature matrix
+ // Get covariance matrix by inverting the curvature matrix. Note that
+ // the curvature matrix is for the unscaled parameters, hence we have
+ // to scale them before we can proceed.
m_covariance = likelihood.curvature()->invert();
- }
+ // Log covariance matrix
+ if (logTerse()) {
+ log << m_covariance << std::endl;
+ }
+
+ } // endif: covariance computation needed
// Write header
if (logTerse()) {
log << std::endl;
- log.header1("Generate butterfly");
+ log.header1("Prepare butterfly computation");
+ }
+
+ // Find parameter indices. We do this by the memory location of the
+ // model parameters.
+ GModelSky* skymodel = dynamic_cast(models[m_srcname]);
+ GModelPar* par_prefactor = &(*skymodel)["Prefactor"];
+ GModelPar* par_index = &(*skymodel)["Index"];
+ int inx_prefactor = -1;
+ int inx_index = -1;
+ for (int i = 0; i < models.npars(); ++i) {
+ if (models.pars()[i] == par_prefactor) {
+ inx_prefactor = i;
+ }
+ if (models.pars()[i] == par_index) {
+ inx_index = i;
+ }
}
- // Initialise dummy time to evaluate spectral model
- GTime time = GTime();
+ // Check indices
- // Initialise vector of gradients
- GVector grad = GVector(m_obs.models().npars());
+ // Store model values and scales
+ double prefactor_mean = (*skymodel)["Prefactor"].value();
+ double index_mean = (*skymodel)["Index"].value();
+ double pivot_mean = (*skymodel)["PivotEnergy"].value();
+ double prefactor_scale = (*skymodel)["Prefactor"].scale();
+ double index_scale = (*skymodel)["Index"].scale();
+ GEnergy pivot(pivot_mean, "MeV");
- // Loop over energy bins positions
- for (int i = 0; i < m_ebounds.size(); ++i) {
+ // Show parameter indices
+ if (logTerse()) {
+ log << gammalib::parformat("Prefactor");
+ log << prefactor_mean << " ph/cm2/s/MeV (";
+ log << inx_prefactor << ")" << std::endl;
+ log << gammalib::parformat("Index");
+ log << index_mean << " (" << inx_index << ")" << std::endl;
+ log << gammalib::parformat("Pivot energy");
+ log << pivot << std::endl;
+ }
- // Initialise number of gradients
- int num_gradient = 0;
+ // Extract covariance elements
+ double cov_pp = m_covariance(inx_prefactor,inx_prefactor);
+ double cov_pg = m_covariance(inx_prefactor,inx_index);
+ double cov_gg = m_covariance(inx_index,inx_index);
- // Get the energy of current bin
- GEnergy energy = m_ebounds.elogmean(i);
+ // Initialise min and max intensities
+ m_min_intensities.assign(m_ebounds.size(), 1e30);
+ m_max_intensities.assign(m_ebounds.size(), 0.0);
- // Header for verbose logging
- if (logExplicit()) {
- std::string msg = "Computing butterfly for bin number "+
- gammalib::str(i)+" at "+energy.print();
- log << std::endl;
- log.header2(msg);
- }
+ // Compute eigenvectors and eigenvalues
+ double lambda1;
+ double lambda2;
+ GVector vector1(2);
+ GVector vector2(2);
+ eigenvectors(cov_pp, cov_pg, cov_pg, cov_gg,
+ &lambda1, &lambda2, &vector1, &vector2);
- // Initialise flux value
- double model_flux = 0.0;
+ // Show eigenvectors and eigenvalues
+ if (logTerse()) {
+ log << gammalib::parformat("Eigenvalue 1") << lambda1 << std::endl;
+ log << gammalib::parformat("Eigenvalue 2") << lambda2 << std::endl;
+ log << gammalib::parformat("Eigenvector 1") << vector1 << std::endl;
+ log << gammalib::parformat("Eigenvector 2") << vector2 << std::endl;
+ }
- // Loop over models
- for (int j = 0; j < models.size(); ++j) {
+ // Confidence scaling
+ double sigma = gammalib::erfinv(m_confidence) * gammalib::sqrt_two;
+ double scale = (sigma*sigma);
- // If source name is found, retrieve the nominal flux
- if (models[j]->name() == m_srcname) {
+ // Show confidence stuff
+ if (logTerse()) {
+ log << gammalib::parformat("Confidence level");
+ log << m_confidence << std::endl;
+ log << gammalib::parformat("Corresponding scaling");
+ log << scale << std::endl;
+ }
- // 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);
+ // Compute minor and major axes of error ellipse
+ double major = scale*std::sqrt(lambda1);
+ double minor = scale*std::sqrt(lambda2);
- }
+ // Write header
+ if (logTerse()) {
+ log << std::endl;
+ log.header1("Generate butterfly");
+ }
- // Loop over model parameters
- for (int i = 0; i < models[j]->size(); ++i) {
+ // Walk around the error ellipse
+ int nt = 360;
+ double dt = gammalib::twopi/double(nt);
+ double t = 0.0;
+ for (int i = 0; i < nt; ++i, t += dt) {
+
+ // Pre-compute sin and cosine
+ double sin_t = std::sin(t);
+ double cos_t = std::cos(t);
+
+ // Compute prefactor and index
+ double prefactor = (major * cos_t * vector1[0] +
+ minor * sin_t * vector2[0]) * prefactor_scale +
+ prefactor_mean;
+ double index = (major * cos_t * vector1[1] +
+ minor * sin_t * vector2[1]) * index_scale +
+ index_mean;
+
+ // Show prefactor and index
+ if (logTerse()) {
+ log << gammalib::parformat("Angle "+gammalib::str(t*gammalib::rad2deg));
+ log << prefactor << " ph/cm2/s/MeV; ";
+ log << index << std::endl;
+ }
- // Get model parameter
- GModelPar parameter = models[j]->at(i);
+ // Setup power law model
+ GModelSpectralPlaw plaw(prefactor, index, pivot);
- // Set gradient value to the vector
- grad[num_gradient] = parameter.gradient();
+ // Evaluate power law at each energy and extract minimum and
+ // maximum intensity
+ for (int k = 0; k < m_ebounds.size(); ++k) {
- // increment vector index
- num_gradient++;
+ // Get the energy of current bin
+ GTime time;
+ GEnergy energy = m_ebounds.elogmean(k);
- } // endfor: Looped over model parameters
+ // Evaluate power law
+ double intensity = plaw.eval(energy, time);
+
+ // Get extremes
+ if (intensity < m_min_intensities[k]) {
+ m_min_intensities[k] = intensity;
+ }
+ if (intensity > m_max_intensities[k]) {
+ m_max_intensities[k] = intensity;
+ }
+
+ } // endfor: looped over all energies
+
+ } // endfor: looped over all angles
+
+ // Write header
+ if (logTerse()) {
+ log << std::endl;
+ log.header1("Results");
+ }
- } // endfor: Looped over models
+ // Setup power law model
+ GModelSpectralPlaw plaw(prefactor_mean, index_mean, pivot);
- // Multiply covariance to the gradient vector
- GVector vector = m_covariance * grad;
+ // Show results
+ for (int k = 0; k < m_ebounds.size(); ++k) {
- // Get the error from the scalar product
- double error = std::sqrt(grad * vector);
+ // Get the energy of current bin
+ GTime time;
+ GEnergy energy = m_ebounds.elogmean(k);
- // Multiply in scale of spectral parameter
- error *= m_par_scale;
+ // Evaluate power law
+ double intensity = plaw.eval(energy, time);
- // Store flux, value and energy for saving
- m_fluxes.push_back(model_flux);
+ // Store results
m_energies.push_back(energy.MeV());
- m_errors.push_back(error);
+ m_intensities.push_back(intensity);
// Log information
- if (logExplicit()) {
- log << " Flux at ";
- log << energy.print();
- log << ": ";
- log << model_flux;
- log << " +- ";
- log << error;
- log << " ph/cm2/s/MeV";
+ if (logTerse()) {
+ log << gammalib::parformat("Intensity at "+energy.print());
+ log << intensity;
+ log << " (";
+ log << m_min_intensities[k];
+ log << ", ";
+ log << m_max_intensities[k];
+ log << ") ph/cm2/s/MeV";
+ log << std::endl;
}
- } //endfor: loop over energy bin
+ } // endfor
// Restore energy dispersion flag for all CTA observations
for (int i = 0; i < m_obs.size(); ++i) {
@@ -384,14 +507,15 @@ void ctbutterfly::save(void)
// Get output filename
m_outfile = (*this)["outfile"].filename();
- // Create CSV table with 3 columns
- GCsv table(m_fluxes.size(), 3);
+ // Create CSV table with 4 columns
+ GCsv table(m_energies.size(), 4);
// Fill CSV table
- for (int i = 0; i < m_fluxes.size(); ++i) {
+ for (int i = 0; i < m_energies.size(); ++i) {
table.real(i, 0, m_energies[i]);
- table.real(i, 1, m_fluxes[i]);
- table.real(i, 2, m_errors[i]);
+ table.real(i, 1, m_intensities[i]);
+ table.real(i, 2, m_min_intensities[i]);
+ table.real(i, 3, m_max_intensities[i]);
}
// Save CSV table
@@ -415,18 +539,25 @@ void ctbutterfly::init_members(void)
{
// Initialise members
m_srcname.clear();
- m_outfile.clear();
- m_ebounds.clear();
+ m_confidence = 0.68;
+ m_max_iter = 50;
m_apply_edisp = false;
+ m_refit = false;
+ m_ebounds.clear();
+ m_outfile.clear();
// Initialise protected members
- m_skymodel = NULL;
- m_par_scale = 0.0;
+ m_opt.clear();
m_obs.clear();
m_covariance.clear();
m_energies.clear();
- m_fluxes.clear();
- m_errors.clear();
+ m_intensities.clear();
+ m_min_intensities.clear();
+ m_max_intensities.clear();
+
+ // Set optimizer parameters
+ m_opt.max_iter(m_max_iter);
+ m_opt.max_stalls(10);
// Return
return;
@@ -441,20 +572,22 @@ void ctbutterfly::init_members(void)
void ctbutterfly::copy_members(const ctbutterfly& app)
{
// Copy attributes
- m_srcname = app.m_srcname;
- m_outfile = app.m_outfile;
- m_ebounds = app.m_ebounds;
+ m_srcname = app.m_srcname;
+ m_confidence = app.m_confidence;
+ m_max_iter = app.m_max_iter;
m_apply_edisp = app.m_apply_edisp;
-
+ m_refit = app.m_refit;
+ m_ebounds = app.m_ebounds;
+ m_outfile = app.m_outfile;
+
// Copy protected members
- m_obs = app.m_obs;
- m_covariance = app.m_covariance;
- m_energies = app.m_energies;
- m_fluxes = app.m_fluxes;
- m_errors = app.m_errors;
-
- // Extract model parameter
- get_model_parameter();
+ m_obs = app.m_obs;
+ m_opt = app.m_opt;
+ m_covariance = app.m_covariance;
+ m_energies = app.m_energies;
+ m_intensities = app.m_intensities;
+ m_min_intensities = app.m_min_intensities;
+ m_max_intensities = app.m_max_intensities;
// Return
return;
@@ -513,15 +646,9 @@ void ctbutterfly::get_parameters(void)
// Get name of test source and check container for this name
m_srcname = (*this)["srcname"].string();
- if (!m_obs.models().contains(m_srcname)) {
- 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);
- }
- // Get relevant model and parameter for model computation
- get_model_parameter();
+ // Check model name and type
+ check_model();
// Create energy boundaries from user parameters
m_ebounds = create_ebounds();
@@ -539,6 +666,10 @@ void ctbutterfly::get_parameters(void)
// Read energy dispersion flag
m_apply_edisp = (*this)["edisp"].boolean();
+ // Read other parameters
+ m_confidence = (*this)["confidence"].real();
+ m_refit = (*this)["refit"].boolean();
+
// Optionally read ahead parameters so that they get correctly
// dumped into the log file
if (read_ahead()) {
@@ -551,62 +682,114 @@ void ctbutterfly::get_parameters(void)
/***********************************************************************//**
- * @brief Get application parameters
+ * @brief Check if sky model is valid
*
* @exception GException::invalid_value
- * Did not find a valid model parameter
+ * Did not find a valid model
*
- * Extracts a pointer to the sky model (m_skymodel) and a pointer to the
- * relevant model parameter (m_model_par) from the model container.
+ * Checks whether the specified sky model exists and whether it is a power
+ * law.
***************************************************************************/
-void ctbutterfly::get_model_parameter(void)
+void ctbutterfly::check_model(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) {
+ GModels& models = const_cast(m_obs.models());
+ GModelSky* model = dynamic_cast(models[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 "
"butterfly computation.";
- throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
+ throw GException::invalid_value(G_CHECK_MODEL, 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);
+ // Check that the spectral model is a power law
+ if (model->spectral()->type() != "PowerLaw") {
+ std::string msg = "\""+model->spectral()->type()+"\" cannot be "
+ "used as spectral model for an butterfly "
+ "computation. Please specify a power law model.";
+ throw GException::invalid_value(G_CHECK_MODEL, msg);
}
- // Set scale value from parameter
- m_par_scale = modelpar->scale();
-
} // endif: source model existed
+ // Otherwise throw an exception
+ else {
+ 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_CHECK_MODEL, msg);
+ }
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Compute normalized eigenvectors and eigenvalues
+ *
+ * @param[in] a Covariance matrix element (0,0)
+ * @param[in] b Covariance matrix element (0,1)
+ * @param[in] c Covariance matrix element (1,0)
+ * @param[in] d Covariance matrix element (1,1)
+ * @param[out] lambda1 First eigenavlue
+ * @param[out] lambda2 Second eigenavlue
+ * @param[out] vector1 First eigenvector (normalized)
+ * @param[out] vector2 Second eigenvector (normalized)
+ *
+ * Computes the eigenvalues and eigenvectors for a 2x2 matrix. See
+ * http://www.math.harvard.edu/archive/21b_fall_04/exhibits/2dmatrices/index.html
+ ***************************************************************************/
+void ctbutterfly::eigenvectors(const double& a,
+ const double& b,
+ const double& c,
+ const double& d,
+ double* lambda1,
+ double* lambda2,
+ GVector* vector1,
+ GVector* vector2)
+{
+ // Compute trace and determinant
+ double trace = a + d;
+ double det = a*d - b*c;
+
+ // Compute eigenvalues
+ double arg1 = 0.5*trace;
+ double arg2 = std::sqrt(0.25*trace*trace - det);
+ *lambda1 = arg1 + arg2;
+ *lambda2 = arg1 - arg2;
+
+ // Compute eigenvectors
+ double vector1_1;
+ double vector1_2;
+ double vector2_1;
+ double vector2_2;
+ if (c != 0.0) {
+ (*vector1)[0] = *lambda1 - det;
+ (*vector1)[1] = c;
+ (*vector2)[0] = *lambda2 - det;
+ (*vector2)[1] = c;
+ }
+ else if (b != 0.0) {
+ (*vector1)[0] = b;
+ (*vector1)[1] = *lambda1 - a;
+ (*vector2)[0] = b;
+ (*vector2)[1] = *lambda2 - a;
+ }
+ else {
+ (*vector1)[0] = 1.0;
+ (*vector1)[1] = 0.0;
+ (*vector2)[0] = 0.0;
+ (*vector2)[1] = 1.0;
+ }
+
+ // Normalize eigenvectors
+ *vector1 /= norm(*vector1);
+ *vector2 /= norm(*vector2);
+
// Return
return;
}
diff --git a/src/ctbutterfly/ctbutterfly.hpp b/src/ctbutterfly/ctbutterfly.hpp
index c23ba12..e1eb0af 100644
--- a/src/ctbutterfly/ctbutterfly.hpp
+++ b/src/ctbutterfly/ctbutterfly.hpp
@@ -78,22 +78,33 @@ protected:
void copy_members(const ctbutterfly& app);
void free_members(void);
void get_parameters(void);
- void get_model_parameter(void);
+ void check_model(void);
+ void eigenvectors(const double& a,
+ const double& b,
+ const double& c,
+ const double& d,
+ double* lambda1,
+ double* lambda2,
+ GVector* vector1,
+ GVector* vector2);
// User parameters
- std::string m_srcname; //!< Name of source to compute butterfly
- std::string m_outfile; //!< Output ascii file
- GEbounds m_ebounds; //!< Energy binning definition
- bool m_apply_edisp; //!< Apply energy dispersion?
+ std::string m_srcname; //!< Name of source to compute butterfly
+ double m_confidence; //!< Confidence level
+ int m_max_iter; //!< Maximum number of iterations
+ bool m_apply_edisp; //!< Apply energy dispersion?
+ bool m_refit; //!< Do refit?
+ GEbounds m_ebounds; //!< Energy binning definition
+ std::string m_outfile; //!< Output ascii file
// 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
+ GOptimizerLM m_opt; //!< Optimizer
+ GObservations m_obs; //!< Observation container
+ GMatrixSparse m_covariance; //!< Covariance matrix
+ std::vector m_energies; //!< Energy values
+ std::vector m_intensities; //!< Power law intensity
+ std::vector m_min_intensities; //!< Minimum intensities
+ std::vector m_max_intensities; //!< Maximum intensities
};
/***********************************************************************//**
diff --git a/src/ctbutterfly/ctbutterfly.par b/src/ctbutterfly/ctbutterfly.par
index aa62f73..7af995b 100644
--- a/src/ctbutterfly/ctbutterfly.par
+++ b/src/ctbutterfly/ctbutterfly.par
@@ -17,17 +17,19 @@
#
# General parameters
#===================
-inobs, f, a, "events.fits",,,"Input event list, counts cube or observation definition file"
-inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
-srcname, s, a, "Crab",,,"Source of interest"
-expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
-psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
-bkgcube, f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
-caldb, s, a, "prod2",,, "Calibration database"
-irf, s, a, "South_50h",,, "Instrument response function"
-edisp, b, h, no,,, "Apply energy dispersion?"
-outfile, f, a, "butterfly.txt",,, "Output ascii file"
-matrix, f, h, "NONE",,, "Input covariance Matrix FITS file"
+inobs, f, a, "events.fits",,,"Input event list, counts cube or observation definition file"
+inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
+srcname, s, a, "Crab",,,"Source of interest"
+expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
+psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
+bkgcube, f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
+caldb, s, a, "prod2",,, "Calibration database"
+irf, s, a, "South_50h",,, "Instrument response function"
+edisp, b, h, no,,, "Apply energy dispersion?"
+outfile, f, a, "butterfly.txt",,, "Output ascii file"
+refit, b, h, no,,, "Do refitting?"
+confidence, r, h, 0.68,0.0,1.0, "Confidence level"
+matrix, f, h, "NONE",,, "Input covariance Matrix FITS file"
#
# Energy binning parameters
diff --git a/src/cterror/cterror.par b/src/cterror/cterror.par
index ea6e626..8def721 100644
--- a/src/cterror/cterror.par
+++ b/src/cterror/cterror.par
@@ -6,7 +6,7 @@
# ------------------------------------------------------------------------- #
# #
# File: cterror.par #
-# Version: 01-00-00 #
+# Version: 1.0.0 #
# #
# Author: Florent Forest #
# florent.forest@isae.fr #