Feature #3386

Implement GModel for Spectral Line

Added by Hatlen Eirik over 3 years ago. Updated almost 2 years ago.

Status:ClosedStart date:10/14/2020
Priority:NormalDue date:
Assigned To:Hatlen Eirik % Done:

100%

Category:-
Target version:2.0.0
Duration:

Description

Suggestion to implement a “Dirac-Delta” spectral model.

As a spectral line is a “smoking gun” signature from dark matter estimating the resulting convoluted spectrum is essential to any analysis looking for dark matter lines.

line-ctools.pdf (145 KB) Preview Hatlen Eirik , 04/22/2021 10:55 PM


Recurrence

No recurrence.

History

#1 Updated by Knödlseder Jürgen about 3 years ago

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.

#2 Updated by Hatlen Eirik almost 3 years ago

  • % Done changed from 0 to 90

#3 Updated by Hatlen Eirik almost 3 years ago

  • % Done changed from 90 to 100

#4 Updated by Hatlen Eirik almost 3 years ago

  • % Done changed from 100 to 90

#5 Updated by Hatlen Eirik almost 3 years ago

  • Target version set to 2.0.0

#6 Updated by Knödlseder Jürgen almost 3 years ago

  • Status changed from New to In Progress

Hi Eirik, is the feature ready to be merged into the code? If yes, please set the status to “Pull request” and indicate in the body of the thread the name of the branch to me merged in. I will then review the code and take care of the merging.

#7 Updated by Hatlen Eirik almost 3 years ago

3386-Implement-GModelSpectral-line
A special method to GResponse::convolve() that handles the energy boundaries for the special case when modelling a narrow Gaussian (width << eres) with energy resolution.

The method returns one boundary used for the whole energy interval.

Attached is a comparison of ctmodel with a narrow Gaussian with energy dispersion enabled.

#8 Updated by Knödlseder Jürgen almost 3 years ago

  • Status changed from In Progress to Pull request

This looks great, thanks a lot, I will proceed with the merging of the code.

#9 Updated by Knödlseder Jürgen almost 3 years ago

  • Status changed from Pull request to Feedback

I checked the code that you committed. I recognised that all the vector-wise computation got removed in your branch, hence I put the old code back and combined it with your modification. I merged the code into devel. Could you please check that things are working on your side?

#10 Updated by Knödlseder Jürgen almost 2 years ago

  • Status changed from Feedback to Closed
  • % Done changed from 90 to 100

No feedback, things seem to work now.

Also available in: Atom PDF