Justine Devin / gammalib

Commit 33d830a20d14c8d34ab1301d6a5e5ad053dfd3ea

Authored by Jürgen Knödlseder
1 parent 51ba4b63

Use GModelSpatialDiffuseCube::set_mc_cone() method to limit the Monte Carlo samp… 

…ling to the simulation cone.

This assures that the simulated spectrum corresponds indeed to the selected zone in the diffuse map cube, and also considerably speeds up computations for large maps.
Inline Side-by-side
Showing 1 changed file with 11 additions and 16 deletions
src/model/GModelSky.cpp
  View file @33d830a
... ... @@ -839,29 +839,24 @@ GPhotons GModelSky::mc(const double& area,
839 839 GModelSpatialDiffuseCube* cube = dynamic_cast<GModelSpatialDiffuseCube*>(m_spatial);
840 840 if (cube != NULL) {
841 841  
842   - // Allocate node function
843   - GModelSpectralNodes* nodes = new GModelSpectralNodes(cube->spectrum());
  842 + // Set MC cone
  843 + cube->set_mc_cone(dir, radius);
844 844  
  845 + // Allocate node function to replace the spectral component
  846 + GModelSpectralNodes* nodes = new GModelSpectralNodes(cube->spectrum());
  847 + for (int i = 0; i < nodes->nodes(); ++i) {
  848 + GEnergy energy = nodes->energy(i);
  849 + double intensity = nodes->intensity(i);
  850 + double norm = m_spectral->eval(energy, tmin);
  851 + nodes->intensity(i, norm*intensity);
  852 + }
  853 +
845 854 // Signal that node function needs to be de-allocated later
846 855 free_spectral = true;
847 856  
848 857 // Set the spectral model pointer to the node function
849 858 spectral = nodes;
850 859  
851   - // Compute the node function intensities as the product of
852   - // the original intensities multiplied by the spectral model
853   - // evaluated at the energies of the node function.
854   - //TODO The spectral model may be time dependent, but at this
855   - // stage we do not know the times; we thus take the minimum
856   - // time, assuming that the spectral model will not change with
857   - // time.
858   - for (int i = 0; i < nodes->nodes(); ++i) {
859   - GEnergy energy = nodes->energy(i);
860   - double intensity = nodes->intensity(i) *
861   - m_spectral->eval(energy, tmin);
862   - nodes->intensity(i, intensity);
863   - }
864   -
865 860 } // endif: spatial model was a diffuse cube
866 861  
867 862 // Compute flux within [emin, emax] in model from spectral
... ...