diff --git a/inst/cta/src/GCTAEdisp2D.cpp b/inst/cta/src/GCTAEdisp2D.cpp index 7b1eeeb..64a28e2 100644 --- a/inst/cta/src/GCTAEdisp2D.cpp +++ b/inst/cta/src/GCTAEdisp2D.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAEdisp2D.cpp - CTA 2D energy dispersion class * * ----------------------------------------------------------------------- * - * copyright (C) 2015-2016 by Florent Forest * + * copyright (C) 2015-2017 by Florent Forest * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -43,6 +43,8 @@ /* __ Method name definitions ____________________________________________ */ #define G_READ "GCTAEdisp2D::read(GFitsTable&)" +#define G_MC "GCTAEdisp2D::mc(GRan&, double&, double&, double&, double&, "\ + "double&)" /* __ Macros _____________________________________________________________ */ @@ -363,7 +365,7 @@ void GCTAEdisp2D::read(const GFitsTable& table) sums.push_back(sum); } // endfor: looped over Etrue - + } // endfor: looped over theta // Normalize vectors of migration matrix for all etrue and theta. @@ -378,7 +380,7 @@ void GCTAEdisp2D::read(const GFitsTable& table) } } } - + } // endfor: made 2 passes // Set maximum energy dispersion value @@ -500,8 +502,12 @@ void GCTAEdisp2D::save(const GFilename& filename, const bool& clobber) const * @param[in] zenith Zenith angle in Earth system (rad). * @param[in] azimuth Azimuth angle in Earth system (rad). * + * @exception GException::invalid_return_value + * No energy dispersion information found for parameters + * * Draws observed energy value given a true energy @p logEsrc and offset - * angle @p theta. + * angle @p theta. If no energy dispersion information is available the + * method will return the true photon energy. ***************************************************************************/ GEnergy GCTAEdisp2D::mc(GRan& ran, const double& logEsrc, @@ -515,10 +521,24 @@ GEnergy GCTAEdisp2D::mc(GRan& ran, double emin = ebounds.emin().log10TeV(); double emax = ebounds.emax().log10TeV(); + // Throw an exception if no valid energy dispersion information is + // available + if (emin == emax) { + std::string msg = "No valid energy dispersion information available " + "for true photon energy "+ebounds.emin().print()+ + ", offset angle "+ + gammalib::str(theta*gammalib::rad2deg)+" deg " + " and azimuth angle "+ + gammalib::str(phi*gammalib::rad2deg)+" deg. Either " + "provide an energy response matrix covering these " + "parameters or restrict the parameter space."; + throw GException::invalid_return_value(G_MC, msg); + } + // Find energy by rejection method double ewidth = emax - emin; - double logEobs = 0.5*(emin+emax); - if (m_max_edisp > 0.0) { + double logEobs = logEsrc; + if (m_max_edisp > 0.0 && ewidth > 0.0) { double f = 0.0; double ftest = 1.0; while (ftest > f) { diff --git a/inst/cta/src/GCTAResponseIrf.cpp b/inst/cta/src/GCTAResponseIrf.cpp index c23ee91..bc9c5f4 100644 --- a/inst/cta/src/GCTAResponseIrf.cpp +++ b/inst/cta/src/GCTAResponseIrf.cpp @@ -489,11 +489,11 @@ double GCTAResponseIrf::nroi(const GModelSky& model, // If energy dispersion is requested then integrate over the relevant // true photon energies ... if (use_edisp()) { - + // Retrieve true energy boundaries double obsLogEng = obsEng.log10TeV(); GEbounds ebounds = edisp()->ebounds_src(obsLogEng); - + // Loop over all boundaries for (int i = 0; i < ebounds.size(); ++i) { @@ -515,7 +515,7 @@ double GCTAResponseIrf::nroi(const GModelSky& model, emin = std::log(emin); emax = std::log(emax); nroi += integral.romberg(emin, emax); - + } // endif: interval was valid } // endfor: looped over energy intervals @@ -524,7 +524,7 @@ double GCTAResponseIrf::nroi(const GModelSky& model, // ... otherwise evaluate else { - + // No energy dispersion const GEnergy& srcEng = obsEng; @@ -535,7 +535,7 @@ double GCTAResponseIrf::nroi(const GModelSky& model, // Compute response nroi = npred_spatial * npred_spectral * npred_temporal; - + } // endelse: no energy dispersion requested // If required, apply instrument specific model scaling @@ -1152,7 +1152,7 @@ void GCTAResponseIrf::load_aeff(const GFilename& filename) m_aeff = new GCTAAeffArf(filename); } // endif: load as GCTAAeffArf - + else { // Close file @@ -1160,7 +1160,7 @@ void GCTAResponseIrf::load_aeff(const GFilename& filename) // Allocate Aeff from file m_aeff = new GCTAAeff2D(filename); - + } // endelse: load as GCTAAeff2D } // endif: extension name is not empty @@ -1270,13 +1270,13 @@ void GCTAResponseIrf::load_psf(const GFilename& filename) // Check for King profile specific table columns if (table.contains("GAMMA") && table.contains("SIGMA")) { - + // Close FITS file fits.close(); // Allocate King profile PSF m_psf = new GCTAPsfKing(filename); - + } // ... otherwise check for Gaussian profile specific table @@ -1284,36 +1284,36 @@ void GCTAResponseIrf::load_psf(const GFilename& filename) else if (table.contains("SCALE") && table.contains("SIGMA_1") && table.contains("AMPL_2") && table.contains("SIGMA_2") && table.contains("AMPL_3") && table.contains("SIGMA_3")) { - + // Close FITS file fits.close(); // Allocate Gaussian profile PSF m_psf = new GCTAPsf2D(filename); - + } // ... otherwise check for PSF table specific table columns else if (table.contains("RAD_LO") && table.contains("RAD_HI") && table.contains("RPSF")) { - + // Close FITS file fits.close(); // Allocate PSF table m_psf = new GCTAPsfTable(filename); - + } - + // ... otherwise try opening as vector PSF else { - + // Close FITS file fits.close(); // Allocate vector PSF m_psf = new GCTAPsfVector(filename); - + } } // endif: extension name is not empty @@ -1331,7 +1331,7 @@ void GCTAResponseIrf::load_psf(const GFilename& filename) // ... otherwise load file as a performance table else { - + // Allocate a performance table m_psf = new GCTAPsfPerfTable(filename); @@ -1405,24 +1405,24 @@ void GCTAResponseIrf::load_edisp(const GFilename& filename) // Check for 2D migration matrix if (table.contains("MIGRA_LO") && table.contains("MIGRA_HI")) { - + // Close FITS file fits.close(); // Allocate 2D migration matrix m_edisp = new GCTAEdisp2D(filename); - + } // ... otherwise allocate RMF else { - + // Close FITS file fits.close(); // Allocate Gaussian profile PSF m_edisp = new GCTAEdispRmf(filename); - + } } // endif: extension name is not empty @@ -2232,7 +2232,7 @@ double GCTAResponseIrf::irf_ptsrc(const GEvent& event, // Set Photon GPhoton photon(src->dir(), source.energy(), source.time()); - + // Compute IRF double irf = this->irf(event, photon, obs);