Dmitry Malyshev / gammalib

Commit 7249d75bf84d4975887330920fa0b0fd37430d97

Authored by Jürgen Knödlseder
1 parent d52edf8f

Make GModelSky::mc() method save against invalid energy ranges or model cubes th… 

…at are not covered by the RoI.
Inline Side-by-side
Showing 1 changed file with 84 additions and 59 deletions
src/model/GModelSky.cpp
  View file @7249d75
1 1 /***************************************************************************
2 2 * GModelSky.cpp - Sky model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -863,14 +863,16 @@ GPhotons GModelSky::mc(const double& area,
863 863 // function spectral model that is the product of the diffuse
864 864 // cube node function and the spectral model evaluated at the
865 865 // energies of the node function
866   - GModelSpatialDiffuseCube* cube = dynamic_cast<GModelSpatialDiffuseCube*>(m_spatial);
  866 + GModelSpatialDiffuseCube* cube =
  867 + dynamic_cast<GModelSpatialDiffuseCube*>(m_spatial);
867 868 if (cube != NULL) {
868 869  
869 870 // Set MC cone
870 871 cube->set_mc_cone(dir, radius);
871 872  
872 873 // Allocate node function to replace the spectral component
873   - GModelSpectralNodes* nodes = new GModelSpectralNodes(cube->spectrum());
  874 + GModelSpectralNodes* nodes =
  875 + new GModelSpectralNodes(cube->spectrum());
874 876 for (int i = 0; i < nodes->nodes(); ++i) {
875 877 GEnergy energy = nodes->energy(i);
876 878 double intensity = nodes->intensity(i);
... ... @@ -884,11 +886,18 @@ GPhotons GModelSky::mc(const double&amp; area,
884 886 // Set the spectral model pointer to the node function
885 887 spectral = nodes;
886 888  
  889 + // Kluge: if there are no nodes then the spectral->flux method
  890 + // will throw an exception. We therefore set here the use_model
  891 + // flag to false in case that there are no spectral nodes
  892 + if (nodes->nodes() == 0) {
  893 + use_model = false;
  894 + }
  895 +
887 896 } // endif: spatial model was a diffuse cube
888 897  
889 898 // Compute flux within [emin, emax] in model from spectral
890 899 // component (units: ph/cm2/s)
891   - double flux = spectral->flux(emin, emax);
  900 + double flux = (use_model) ? spectral->flux(emin, emax) : 0.0;
892 901  
893 902 // Derive expecting counting rate within simulation surface
894 903 // (units: ph/s)
... ... @@ -902,72 +911,88 @@ GPhotons GModelSky::mc(const double&amp; area,
902 911 std::cout << "norm=" << norm << ")" << std::endl;
903 912 #endif
904 913  
905   - // Get photon arrival times from temporal model
906   - GTimes times = m_temporal->mc(rate, tmin, tmax, ran);
907   -
908   - // Debug option: dump number of times
909   - #if defined(G_DUMP_MC_DETAIL)
910   - std::cout << " Times=" << times.size() << std::endl;
911   - #endif
912   -
913   - // Reserve space for photons
914   - if (times.size() > 0) {
915   - photons.reserve(times.size());
916   - }
917   -
918   - // Loop over photons
919   - for (int i = 0; i < times.size(); ++i) {
920   -
921   - // Debug option: dump photon index
922   - #if defined(G_DUMP_MC_DETAIL)
923   - std::cout << " Photon=" << i << std::endl;
924   - #endif
925   -
926   - // Allocate photon
927   - GPhoton photon;
  914 + // Continue only if rate is positive
  915 + if (rate > 0.0) {
928 916  
929   - // Set photon arrival time
930   - photon.time(times[i]);
  917 + // Get photon arrival times from temporal model
  918 + GTimes times = m_temporal->mc(rate, tmin, tmax, ran);
931 919  
932   - // Debug option: dump time
  920 + // Debug option: dump number of times
933 921 #if defined(G_DUMP_MC_DETAIL)
934   - std::cout << " Time=" << times[i] << std::endl;
  922 + std::cout << " Times=" << times.size() << std::endl;
935 923 #endif
936 924  
937   - // Set photon energy
938   - photon.energy(spectral->mc(emin, emax, photon.time(), ran));
939   -
940   - // Debug option: dump energy
941   - #if defined(G_DUMP_MC_DETAIL)
942   - std::cout << " Energy=" << photon.energy() << std::endl;
943   - #endif
944   -
945   - // Set incident photon direction. If an invalid_return_value
946   - // exception occurs the sky direction returned by the spatial
947   - // Monte Carlo method is invalid and the photon is skipped.
948   - try {
949   - photon.dir(m_spatial->mc(photon.energy(), photon.time(), ran));
950   - }
951   - catch (GException::invalid_return_value) {
952   - continue;
  925 + // Reserve space for photons
  926 + if (times.size() > 0) {
  927 + photons.reserve(times.size());
953 928 }
954 929  
955   - // Debug option: dump direction
956   - #if defined(G_DUMP_MC_DETAIL)
957   - std::cout << " Direction=" << photon.dir() << std::endl;
958   - #endif
959   -
960   - // Append photon
961   - if (dir.dist_deg(photon.dir()) <= radius) {
962   - photons.append(photon);
963   - }
964   -
965   - } // endfor: looped over photons
  930 + // Loop over photons
  931 + for (int i = 0; i < times.size(); ++i) {
  932 +
  933 + // Debug option: dump photon index
  934 + #if defined(G_DUMP_MC_DETAIL)
  935 + std::cout << " Photon=" << i << std::endl;
  936 + #endif
  937 +
  938 + // Allocate photon
  939 + GPhoton photon;
  940 +
  941 + // Set photon arrival time
  942 + photon.time(times[i]);
  943 +
  944 + // Debug option: dump time
  945 + #if defined(G_DUMP_MC_DETAIL)
  946 + std::cout << " Time=" << times[i] << std::endl;
  947 + #endif
  948 +
  949 + // Set photon energy. If an invalid_return_value exception
  950 + // occurs the energy returned by the spectral Monte Carlo
  951 + // method is invalid and the photon is skipped.
  952 + try {
  953 + photon.energy(spectral->mc(emin, emax, photon.time(),
  954 + ran));
  955 + }
  956 + catch (GException::invalid_return_value) {
  957 + continue;
  958 + }
  959 +
  960 + // Debug option: dump energy
  961 + #if defined(G_DUMP_MC_DETAIL)
  962 + std::cout << " Energy=" << photon.energy() << std::endl;
  963 + #endif
  964 +
  965 + // Set incident photon direction. If an invalid_return_value
  966 + // exception occurs the sky direction returned by the
  967 + // spatial Monte Carlo method is invalid and the photon is
  968 + // skipped.
  969 + try {
  970 + photon.dir(m_spatial->mc(photon.energy(), photon.time(),
  971 + ran));
  972 + }
  973 + catch (GException::invalid_return_value) {
  974 + continue;
  975 + }
  976 +
  977 + // Debug option: dump direction
  978 + #if defined(G_DUMP_MC_DETAIL)
  979 + std::cout << " Direction=" << photon.dir() << std::endl;
  980 + #endif
  981 +
  982 + // Append photon
  983 + if (dir.dist_deg(photon.dir()) <= radius) {
  984 + photons.append(photon);
  985 + }
  986 +
  987 + } // endfor: looped over photons
  988 +
  989 + } // endif: rate was positive
966 990  
967 991 // Free spectral model if required
968 992 if (free_spectral) delete spectral;
969 993  
970 994 } // endif: model was used
  995 +
971 996 } // endif: model was valid
972 997  
973 998 // Return photon list
... ... @@ -978,7 +1003,7 @@ GPhotons GModelSky::mc(const double&amp; area,
978 1003 /***********************************************************************//**
979 1004 * @brief Print model information
980 1005 *
981   - * @param[in] chatter Chattiness (defaults to NORMAL).
  1006 + * @param[in] chatter Chattiness.
982 1007 * @return String containing model information.
983 1008 ***************************************************************************/
984 1009 std::string GModelSky::print(const GChatter& chatter) const
... ...