GammaLib

You can in principle do this already by using an isotrop diffuse model as background component. But I agree that this will be time consuming as it involves a convolution.

Ah, I understand your point and haven’t thought about this option. But you are right: especially for a large datasets this could cause computing problems. I am also thinking about the binned/stacked analysis where we need to compute a background cube first. Currently, only background models and no sky models are used in ctbkgcube.
Do you think creating GCTABackground2D is a reasonable approach then?

I’m not sure we need a new GCTABackground... class for that, won’t it be sufficient to create a new model, e.g. GCTAModelAeffBackground (which should be pretty similar to GCTAModelIrfBackground)?

Actually, I had this in mind as well (with exactly that classname ). Going through the code however, I realised that the GCTAModelAeffBackground::mc-method would require some thinking as GCTAAeff doesn’t contain a method to create randomised event directions. But I guess this can be easily adapted from GCTABackground?

Mayer Michael wrote:

Actually, I had this in mind as well (with exactly that classname ). Going through the code however, I realised that the GCTAModelAeffBackground::mc-method would require some thinking as GCTAAeff doesn’t contain a method to create randomised event directions. But I guess this can be easily adapted from GCTABackground?

Indeed, the mc method would need to be coded explicitly in GCTAModelAeffBackground.

Ok, I will give it a try ...

And I agree, using GCTAModelAeffBackground allows for more flexibility if the effective area format gets changed or extended eventually.

As we do not need this for the other functionalities I would not do this (it creates an unnecessary burden to implement something for all classes that is generally not needed).

Can’t we copy basically the mc method from GCTABackground3D?

But I see your problem: you can’t have an computation stack as we’re living in model space. But we could implement a simple rejection method, right? It would only require that we know the maximum Aeff value.

Can’t we copy basically the mc method from GCTABackground3D? But I see your problem: you can’t have an computation stack as we’re living in model space. But we could implement a simple rejection method, right? It would only require that we know the maximum Aeff value.

You are right. I was wondering how to sample the Aeff function to find the maximum value. I guess for the energy range we could simply use the GEbounds that come with the observation, while for the offset, we could use the GCTAPointing and the GCTARoi to find the maximum range. The only thing we need is then a hard-coded sampling granularity in energy and offset.

Mayer Michael wrote:

Can’t we copy basically the mc method from GCTABackground3D? But I see your problem: you can’t have an computation stack as we’re living in model space. But we could implement a simple rejection method, right? It would only require that we know the maximum Aeff value.

You are right. I was wondering how to sample the Aeff function to find the maximum value. I guess for the energy range we could simply use the GEbounds that come with the observation, while for the offset, we could use the GCTAPointing and the GCTARoi to find the maximum range. The only thing we need is then a hard-coded sampling granularity in energy and offset.

Why do you need any granularity for a rejection method?

Sorry, I wasn’t expressing myself clearly. I need the granularity for finding the maximum effective area:

 // Initialise maximum effective area
double max_aeff = 0.0;

// Granularity
int m_onodes = 20;

// Loop over offsets to find maximum effective area
for (int j = 0; j < m_onodes; ++j) {
     double theta = j * max_theta / (double(m_onodes) - 1.0);
     double aeff_theta = (*aeff)(logE, theta)
     if (aeff_theta > max_aeff) {
         max_aeff = aeff_theta;
     }
}

For randomly drawing an event energy, I would also use the GModelSpectralNodes::mc method, where I also need to define a number of nodes.

I guess we can use the granularity of the Aeff matrix. Values are interpolated bilinearily, hence we have to simply find the largest value by looping for a given energy over all offaxis angles.

I guess we can use the granularity of the Aeff matrix. Values are interpolated bilinearily, hence we have to simply find the largest value by looping for a given energy over all offaxis angles.

That is true for the case of GCTAAeff2D. However, in the code, we do not know what kind of GCTAAeff we are dealing with (GCTAAeffArf, GCTAAeffPerfTable, GCTAAeff2D, or any future Aeff format). In principle, it should work for any given effective area which implements operator()(logE, theta, phi, zenith, azimuth, etrue), right? Or would you limit this feature to GCTAAeff2D?

This is true. I’m wondering whether we may implement an GCTAAeff::max(GEnergy) method (or something like that) so that we do not have to recompute the value every time.

I was looking into the code and I was wondering whether the Monte Carlo simulations of the background model were in fact correct. I see that you use:

offset = std::sqrt(ran.uniform()) * max_theta;
phi    = ran.uniform() * gammalib::twopi;

cosrad = std::cos(max_theta);
offset = std::acos(1.0 - ran.uniform() * (1.0 - cosrad));
phi    = ran.uniform() * gammalib::twopi;

I overlooked the sqrt-term in your equation, which corresponds to the small angle approximation of my formula. Sorry for that.

I nevertheless propose to use the correct trigonometric formula.

I agree with you regarding the correct trigonometric formula. Thanks for merging.