GammaLib

I think that spectral line fitting only makes sense with energy dispersion enabled.

Limiting to this case it is probably not very difficult to implement. There is already a spectral model GModelSpectralGauss that can be used for this purpose.

What is needed is to add some code in GResponse::convolve() where the true energy interval is determined:

// Retrieve true energy boundaries
GEbounds ebounds = this->ebounds(event.energy());

// Initialise vector array
GVector array(size);

// Loop over all boundaries
for (int i = 0; i < ebounds.size(); ++i) {

    // Get true energy boundaries in MeV
    double etrue_min = ebounds.emin(i).MeV();
    double etrue_max = ebounds.emax(i).MeV();

    // Continue only if valid
    if (etrue_max > etrue_min) {

        // Setup integration function
        edisp_kerns integrand(this, &obs, &model, &event, srcTime, grad);
        GIntegrals  integral(&integrand);

        // Set number of iterations
        integral.fixed_iter(iter);

        // Do Romberg integration
        array += integral.romberg(etrue_min, etrue_max);

    } // endif: interval was valid

} // endfor: looped over intervals

etrue_min and etrue_max is the true energy interval over which the energy dispersion integration is done for a given reconstructed energy.

One would simply need to add some code that computes the validity true energy interval for a spectral model. This should be done in some private method so that the same code can also be used elsewhere. The method could for example be:

GEbounds GResponse::ebounds_model(const GModelSky& model)

For all spectral except of GModelSpectralGauss the method would return an empty energy boundary container. For GModelSpectralGauss the method should return a single energy boundary that is for example 5 sigma around the line energy.

Some code needs then to be added to limit the energy range for the integration:

// Retrieve true energy boundaries
GEbounds ebounds = this->ebounds(event.energy());

// Initialise vector array
GVector array(size);

// Loop over all boundaries
for (int i = 0; i < ebounds.size(); ++i) {

    // Get true energy boundaries in MeV
    double etrue_min = ebounds.emin(i).MeV();
    double etrue_max = ebounds.emax(i).MeV();

    // Get model energy boundaries
    GEbounds ebounds_model = this->ebounds_model(model);

    // Limit energy boundaries
    if (ebounds_model.size() > 0) {
        if (ebounds_model.emin(i).MeV() > etrue_min) {
            etrue_min = ebounds_model.emin(i).MeV();
        }
        if (ebounds_model.emax(i).MeV() < etrue_max) {
            etrue_max = ebounds_model.emax(i).MeV();
        }
    }

    // Continue only if valid
    if (etrue_max > etrue_min) {

        ...

    } // endif: interval was valid

} // endfor: looped over intervals

Note that the call to GResponse::ebounds_model should be done outside the event loop to save computing time. It was only put for illustration within the event loop.