Dmitry Malyshev / gammalib

Commit d8c52f9d4aed2104696eeaa912f05eea3eddf487

Authored by Jürgen Knödlseder
1 parent c97a75db

Throw GException::invalid_return_value exception in GCTAEdisp2D is no energy dis… 

…persion information is available
Inline Side-by-side
Showing 2 changed files with 48 additions and 28 deletions
inst/cta/src/GCTAEdisp2D.cpp
  View file @d8c52f9
1 1 /***************************************************************************
2 2 * GCTAEdisp2D.cpp - CTA 2D energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2015-2016 by Florent Forest *
  4 + * copyright (C) 2015-2017 by Florent Forest *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -43,6 +43,8 @@
43 43  
44 44 /* __ Method name definitions ____________________________________________ */
45 45 #define G_READ "GCTAEdisp2D::read(GFitsTable&)"
  46 +#define G_MC "GCTAEdisp2D::mc(GRan&, double&, double&, double&, double&, "\
  47 + "double&)"
46 48  
47 49 /* __ Macros _____________________________________________________________ */
48 50  
... ... @@ -363,7 +365,7 @@ void GCTAEdisp2D::read(const GFitsTable& table)
363 365 sums.push_back(sum);
364 366  
365 367 } // endfor: looped over Etrue
366   -
  368 +
367 369 } // endfor: looped over theta
368 370  
369 371 // Normalize vectors of migration matrix for all etrue and theta.
... ... @@ -378,7 +380,7 @@ void GCTAEdisp2D::read(const GFitsTable& table)
378 380 }
379 381 }
380 382 }
381   -
  383 +
382 384 } // endfor: made 2 passes
383 385  
384 386 // Set maximum energy dispersion value
... ... @@ -500,8 +502,12 @@ void GCTAEdisp2D::save(const GFilename& filename, const bool& clobber) const
500 502 * @param[in] zenith Zenith angle in Earth system (rad).
501 503 * @param[in] azimuth Azimuth angle in Earth system (rad).
502 504 *
  505 + * @exception GException::invalid_return_value
  506 + * No energy dispersion information found for parameters
  507 + *
503 508 * Draws observed energy value given a true energy @p logEsrc and offset
504   - * angle @p theta.
  509 + * angle @p theta. If no energy dispersion information is available the
  510 + * method will return the true photon energy.
505 511 ***************************************************************************/
506 512 GEnergy GCTAEdisp2D::mc(GRan& ran,
507 513 const double& logEsrc,
... ... @@ -515,10 +521,24 @@ GEnergy GCTAEdisp2D::mc(GRan& ran,
515 521 double emin = ebounds.emin().log10TeV();
516 522 double emax = ebounds.emax().log10TeV();
517 523  
  524 + // Throw an exception if no valid energy dispersion information is
  525 + // available
  526 + if (emin == emax) {
  527 + std::string msg = "No valid energy dispersion information available "
  528 + "for true photon energy "+ebounds.emin().print()+
  529 + ", offset angle "+
  530 + gammalib::str(theta*gammalib::rad2deg)+" deg "
  531 + " and azimuth angle "+
  532 + gammalib::str(phi*gammalib::rad2deg)+" deg. Either "
  533 + "provide an energy response matrix covering these "
  534 + "parameters or restrict the parameter space.";
  535 + throw GException::invalid_return_value(G_MC, msg);
  536 + }
  537 +
518 538 // Find energy by rejection method
519 539 double ewidth = emax - emin;
520   - double logEobs = 0.5*(emin+emax);
521   - if (m_max_edisp > 0.0) {
  540 + double logEobs = logEsrc;
  541 + if (m_max_edisp > 0.0 && ewidth > 0.0) {
522 542 double f = 0.0;
523 543 double ftest = 1.0;
524 544 while (ftest > f) {
... ...
inst/cta/src/GCTAResponseIrf.cpp
  View file @d8c52f9
... ... @@ -489,11 +489,11 @@ double GCTAResponseIrf::nroi(const GModelSky& model,
489 489 // If energy dispersion is requested then integrate over the relevant
490 490 // true photon energies ...
491 491 if (use_edisp()) {
492   -
  492 +
493 493 // Retrieve true energy boundaries
494 494 double obsLogEng = obsEng.log10TeV();
495 495 GEbounds ebounds = edisp()->ebounds_src(obsLogEng);
496   -
  496 +
497 497 // Loop over all boundaries
498 498 for (int i = 0; i < ebounds.size(); ++i) {
499 499  
... ... @@ -515,7 +515,7 @@ double GCTAResponseIrf::nroi(const GModelSky&amp; model,
515 515 emin = std::log(emin);
516 516 emax = std::log(emax);
517 517 nroi += integral.romberg(emin, emax);
518   -
  518 +
519 519 } // endif: interval was valid
520 520  
521 521 } // endfor: looped over energy intervals
... ... @@ -524,7 +524,7 @@ double GCTAResponseIrf::nroi(const GModelSky&amp; model,
524 524  
525 525 // ... otherwise evaluate
526 526 else {
527   -
  527 +
528 528 // No energy dispersion
529 529 const GEnergy& srcEng = obsEng;
530 530  
... ... @@ -535,7 +535,7 @@ double GCTAResponseIrf::nroi(const GModelSky&amp; model,
535 535  
536 536 // Compute response
537 537 nroi = npred_spatial * npred_spectral * npred_temporal;
538   -
  538 +
539 539 } // endelse: no energy dispersion requested
540 540  
541 541 // If required, apply instrument specific model scaling
... ... @@ -1152,7 +1152,7 @@ void GCTAResponseIrf::load_aeff(const GFilename&amp; filename)
1152 1152 m_aeff = new GCTAAeffArf(filename);
1153 1153  
1154 1154 } // endif: load as GCTAAeffArf
1155   -
  1155 +
1156 1156 else {
1157 1157  
1158 1158 // Close file
... ... @@ -1160,7 +1160,7 @@ void GCTAResponseIrf::load_aeff(const GFilename&amp; filename)
1160 1160  
1161 1161 // Allocate Aeff from file
1162 1162 m_aeff = new GCTAAeff2D(filename);
1163   -
  1163 +
1164 1164 } // endelse: load as GCTAAeff2D
1165 1165  
1166 1166 } // endif: extension name is not empty
... ... @@ -1270,13 +1270,13 @@ void GCTAResponseIrf::load_psf(const GFilename&amp; filename)
1270 1270  
1271 1271 // Check for King profile specific table columns
1272 1272 if (table.contains("GAMMA") && table.contains("SIGMA")) {
1273   -
  1273 +
1274 1274 // Close FITS file
1275 1275 fits.close();
1276 1276  
1277 1277 // Allocate King profile PSF
1278 1278 m_psf = new GCTAPsfKing(filename);
1279   -
  1279 +
1280 1280 }
1281 1281  
1282 1282 // ... otherwise check for Gaussian profile specific table
... ... @@ -1284,36 +1284,36 @@ void GCTAResponseIrf::load_psf(const GFilename&amp; filename)
1284 1284 else if (table.contains("SCALE") && table.contains("SIGMA_1") &&
1285 1285 table.contains("AMPL_2") && table.contains("SIGMA_2") &&
1286 1286 table.contains("AMPL_3") && table.contains("SIGMA_3")) {
1287   -
  1287 +
1288 1288 // Close FITS file
1289 1289 fits.close();
1290 1290  
1291 1291 // Allocate Gaussian profile PSF
1292 1292 m_psf = new GCTAPsf2D(filename);
1293   -
  1293 +
1294 1294 }
1295 1295  
1296 1296 // ... otherwise check for PSF table specific table columns
1297 1297 else if (table.contains("RAD_LO") && table.contains("RAD_HI") &&
1298 1298 table.contains("RPSF")) {
1299   -
  1299 +
1300 1300 // Close FITS file
1301 1301 fits.close();
1302 1302  
1303 1303 // Allocate PSF table
1304 1304 m_psf = new GCTAPsfTable(filename);
1305   -
  1305 +
1306 1306 }
1307   -
  1307 +
1308 1308 // ... otherwise try opening as vector PSF
1309 1309 else {
1310   -
  1310 +
1311 1311 // Close FITS file
1312 1312 fits.close();
1313 1313  
1314 1314 // Allocate vector PSF
1315 1315 m_psf = new GCTAPsfVector(filename);
1316   -
  1316 +
1317 1317 }
1318 1318  
1319 1319 } // endif: extension name is not empty
... ... @@ -1331,7 +1331,7 @@ void GCTAResponseIrf::load_psf(const GFilename&amp; filename)
1331 1331  
1332 1332 // ... otherwise load file as a performance table
1333 1333 else {
1334   -
  1334 +
1335 1335 // Allocate a performance table
1336 1336 m_psf = new GCTAPsfPerfTable(filename);
1337 1337  
... ... @@ -1405,24 +1405,24 @@ void GCTAResponseIrf::load_edisp(const GFilename&amp; filename)
1405 1405  
1406 1406 // Check for 2D migration matrix
1407 1407 if (table.contains("MIGRA_LO") && table.contains("MIGRA_HI")) {
1408   -
  1408 +
1409 1409 // Close FITS file
1410 1410 fits.close();
1411 1411  
1412 1412 // Allocate 2D migration matrix
1413 1413 m_edisp = new GCTAEdisp2D(filename);
1414   -
  1414 +
1415 1415 }
1416 1416  
1417 1417 // ... otherwise allocate RMF
1418 1418 else {
1419   -
  1419 +
1420 1420 // Close FITS file
1421 1421 fits.close();
1422 1422  
1423 1423 // Allocate Gaussian profile PSF
1424 1424 m_edisp = new GCTAEdispRmf(filename);
1425   -
  1425 +
1426 1426 }
1427 1427  
1428 1428 } // endif: extension name is not empty
... ... @@ -2232,7 +2232,7 @@ double GCTAResponseIrf::irf_ptsrc(const GEvent&amp; event,
2232 2232  
2233 2233 // Set Photon
2234 2234 GPhoton photon(src->dir(), source.energy(), source.time());
2235   -
  2235 +
2236 2236 // Compute IRF
2237 2237 double irf = this->irf(event, photon, obs);
2238 2238  
... ...