@robertazanin - I continued a bit with this and pushed two commits here.
I also remembered that Fabio and Lea started the same exercise during the Paris workshop:
There's the fundamental question here if we should use Sherpa or just continue implementing this ourselves in Gammapy just using Numpy. My feeling is that the Sherpa 3D stuff is very complex, and continuing with the simple 3D model and npred implementation just using Numpy and Gammapy we have started here could be useful (and once available we can discuss which version we want to keep and maintain in Gammapy, either the Sherpa or homegrown or both). So I'd say we continue with this example a bit, and try to split out clearly useful parts into the Gammapy library.
The next steps are to apply the PSF convolution and the EDISP.
@cdeil +1 to reimplement this with numpy and split out useful parts into Gammapy. I think that's what we've also been doing for images and spectra.
PSF convolution works via (pseudo code...):
Applying energy dispersion doesn't exist...