Justine Devin / gammalib

Commit 994f34621c742701d2a94d2918f40f88ab63f220

Authored by Jürgen Knödlseder
1 parent 5873fc2f

Adapt GCTAResponseIrf::load_xxx() methods to new handling of extension names.

Inline Side-by-side
Showing 1 changed file with 224 additions and 94 deletions
inst/cta/src/GCTAResponseIrf.cpp
  View file @994f346
... ... @@ -34,6 +34,7 @@
34 34 #include "GFits.hpp"
35 35 #include "GTools.hpp"
36 36 #include "GMath.hpp"
  37 +#include "GFilename.hpp"
37 38 #include "GIntegral.hpp"
38 39 #include "GCaldb.hpp"
39 40 #include "GSource.hpp"
... ... @@ -1071,17 +1072,21 @@ void GCTAResponseIrf::load(const std::string& rspname)
1071 1072 * @param[in] filename Effective area filename.
1072 1073 *
1073 1074 * @exception GException::file_error
1074   - * File not found.
  1075 + * File or extension not found.
  1076 + *
  1077 + * Loads the effective area from a response file.
1075 1078 *
1076   - * Loads the effective area from a response file. If the file is a FITS file
1077   - * the method first checks whether it contains an "EFFECTIVE AREA" extension.
1078   - * If such an extension is found, a GCTAAeff2D effective area will be
1079   - * allocated. Otherwise, the method checks whether the FITS file contains
1080   - * a "SPECRESP" extension, and if this is the case, the method will allocate
1081   - * a GCTAAeffArf effective area object. In both cases, the file will
1082   - * extract the optional "LO_THRES" and "HI_THRES" save energy thresholds
1083   - * from the FITS file header. If the file is not a FITS file, it will be
1084   - * interpreted as a GCTAAeffPerfTable performance table.
  1079 + * If the file is a FITS file, the method will either use the extension name
  1080 + * specified with the filename, or if no extension name is given, search for
  1081 + * an "EFFECTIVE AREA" or "SPECRESP" extension in the file and open the
  1082 + * corresponding FITS table. If a column named "SPECRESP" is found in the
  1083 + * table, a GCTAAeffArf object will be allocated. Otherwise, a GCTAAeff2D
  1084 + * object will be allocated. In both cases, the method will extract the
  1085 + * optional "LO_THRES" and "HI_THRES" save energy thresholds from the FITS
  1086 + * file header.
  1087 + *
  1088 + * If the file is not a FITS file, it will be interpreted as a
  1089 + * GCTAAeffPerfTable performance table.
1085 1090 ***************************************************************************/
1086 1091 void GCTAResponseIrf::load_aeff(const std::string& filename)
1087 1092 {
... ... @@ -1089,8 +1094,11 @@ void GCTAResponseIrf::load_aeff(const std::string& filename)
1089 1094 if (m_aeff != NULL) delete m_aeff;
1090 1095 m_aeff = NULL;
1091 1096  
  1097 + // Initialise filename
  1098 + GFilename fname(filename);
  1099 +
1092 1100 // Check for existence of file
1093   - if (!gammalib::file_exists_gzip(filename)) {
  1101 + if (!gammalib::file_exists_gzip(fname.filename())) {
1094 1102 std::string msg = "File \""+filename+"\" not found. Please specify "
1095 1103 "a valid effective area response file.";
1096 1104 throw GException::file_error(G_LOAD_AEFF, msg);
... ... @@ -1100,63 +1108,79 @@ void GCTAResponseIrf::load_aeff(const std::string& filename)
1100 1108 if (gammalib::is_fits(filename)) {
1101 1109  
1102 1110 // Open FITS file
1103   - GFits file(filename);
  1111 + GFits file(fname.filename());
  1112 +
  1113 + // Get the extension name. If an extension name has been specified
  1114 + // then use this name, otherwise use either the "EFFECTIVE AREA"
  1115 + // or the "SPECRESP" extension.
  1116 + std::string extname = "";
  1117 + if (fname.has_extname()) {
  1118 + extname = fname.extname();
  1119 + }
  1120 + else {
  1121 + if (file.contains("EFFECTIVE AREA")) {
  1122 + extname = "EFFECTIVE AREA";
  1123 + }
  1124 + else if (file.contains("SPECRESP")) {
  1125 + extname = "SPECRESP";
  1126 + }
  1127 + }
1104 1128  
1105   - // If file contains an "EFFECTIVE AREA" extension then load it
1106   - // as CTA response table
1107   - if (file.contains("EFFECTIVE AREA")) {
  1129 + // Continue only if extension name is not empty
  1130 + if (!extname.empty()) {
1108 1131  
1109   - // Get HDU
1110   - const GFitsHDU* hdu = file.at("EFFECTIVE AREA");
  1132 + // Get FITS table
  1133 + const GFitsTable& table = *file.table(extname);
1111 1134  
1112 1135 // Read save energy thresholds if available
1113   - if (hdu->has_card("LO_THRES")) {
1114   - m_lo_save_thres = hdu->real("LO_THRES");
  1136 + if (table.has_card("LO_THRES")) {
  1137 + m_lo_save_thres = table.real("LO_THRES");
1115 1138 }
1116   - if (hdu->has_card("HI_THRES")) {
1117   - m_hi_save_thres = hdu->real("HI_THRES");
  1139 + if (table.has_card("HI_THRES")) {
  1140 + m_hi_save_thres = table.real("HI_THRES");
1118 1141 }
1119 1142  
1120   - // Close file
1121   - file.close();
1122   -
1123   - // Allocate Aeff from file
1124   - m_aeff = new GCTAAeff2D(filename);
  1143 + // Check for specific table column
  1144 + if (table.contains("SPECRESP")) {
1125 1145  
1126   - }
  1146 + // Close file
  1147 + file.close();
1127 1148  
1128   - // ... else if file contains a "SPECRESP" extension then load it
1129   - // as ARF
1130   - else if (file.contains("SPECRESP")) {
  1149 + // Allocate Aeff from file
  1150 + m_aeff = new GCTAAeffArf(filename);
1131 1151  
1132   - // Get HDU
1133   - const GFitsHDU* hdu = file.at("SPECRESP");
  1152 + } // endif: load as GCTAAeffArf
  1153 +
  1154 + else {
1134 1155  
1135   - // Read save energy thresholds if available
1136   - if (hdu->has_card("LO_THRES")) {
1137   - m_lo_save_thres = hdu->real("LO_THRES");
1138   - }
1139   - if (hdu->has_card("LO_THRES")) {
1140   - m_hi_save_thres = hdu->real("HI_THRES");
1141   - }
  1156 + // Close file
  1157 + file.close();
1142 1158  
1143   - // Close file
1144   - file.close();
  1159 + // Allocate Aeff from file
  1160 + m_aeff = new GCTAAeff2D(filename);
  1161 +
  1162 + } // endelse: load as GCTAAeff2D
1145 1163  
1146   - // Allocate Aeff from file
1147   - m_aeff = new GCTAAeffArf(filename);
  1164 + } // endif: extension name is not empty
1148 1165  
  1166 + // Signal that no extension was found
  1167 + else {
  1168 + std::string msg = "FITS file \""+filename+"\" does not "
  1169 + "contain a valid effective area table. "
  1170 + "Please specify a valid effective area "
  1171 + "response file.";
  1172 + throw GException::file_error(G_LOAD_AEFF, msg);
1149 1173 }
1150 1174  
1151   - }
  1175 + } // endif: file was FITS file
1152 1176  
1153   - // ... otherwise handle ASCII performance tables
  1177 + // ... else try handling an ASCII performance table
1154 1178 else {
1155 1179  
1156 1180 // Allocate a performance table
1157 1181 m_aeff = new GCTAAeffPerfTable(filename);
1158 1182  
1159   - }
  1183 + } // endif: load as GCTAAeffPerfTable
1160 1184  
1161 1185 // Return
1162 1186 return;
... ... @@ -1169,17 +1193,21 @@ void GCTAResponseIrf::load_aeff(const std::string& filename)
1169 1193 * @param[in] filename FITS file name.
1170 1194 *
1171 1195 * @exception GException::file_error
1172   - * File not found.
  1196 + * File or extension not found.
  1197 + *
  1198 + * Loads the point spead function from a response file.
1173 1199 *
1174   - * Loads the point spead function from a response file. If the file is a
1175   - * FITS file, the method first checks whether it contains a
1176   - * "POINT SPREAD FUNCTION" extension. If such an extension is found,
1177   - * either a GCTAPsfKing or a GCTAPsf2D point spread function will be
1178   - * allocated depending on the column names that are found in the extension.
1179   - * If the "POINT SPREAD FUNCTION" extension is not found, the FITS file will
1180   - * be interpreted as a GCTAPsfVector point spread function. If the file is
1181   - * not a FITS file, it will be interpreted as a GCTAPsfPerfTable performance
1182   - * table.
  1200 + * If the file is a FITS file, the method will either use the extension name
  1201 + * specified with the filename, or if no extension name is given, search for
  1202 + * a "POINT SPREAD FUNCTION" or "PSF" extension in the file and open the
  1203 + * corresponding FITS table. If columns named "GAMMA" and "SIGMA" are found
  1204 + * in the table, a GCTAPsfKing object will be allocated. If columns named
  1205 + * "SCALE", "SIGMA_1", "AMPL_2", "SIGMA_2", "AMPL_3" and "SIGMA_3" are found
  1206 + * in the table, a GCTAPsf2D object will be allocated. Otherwise, a
  1207 + * GCTAPsfVector object will be allocated.
  1208 + *
  1209 + * If the file is not a FITS file, it will be interpreted as a
  1210 + * GCTAPsfPerfTable performance table.
1183 1211 ***************************************************************************/
1184 1212 void GCTAResponseIrf::load_psf(const std::string& filename)
1185 1213 {
... ... @@ -1187,8 +1215,11 @@ void GCTAResponseIrf::load_psf(const std::string& filename)
1187 1215 if (m_psf != NULL) delete m_psf;
1188 1216 m_psf = NULL;
1189 1217  
  1218 + // Initialise filename
  1219 + GFilename fname(filename);
  1220 +
1190 1221 // Check for existence of file
1191   - if (!gammalib::file_exists_gzip(filename)) {
  1222 + if (!gammalib::file_exists_gzip(fname.filename())) {
1192 1223 std::string msg = "File \""+filename+"\" not found. Please specify "
1193 1224 "a valid point spread function response file.";
1194 1225 throw GException::file_error(G_LOAD_PSF, msg);
... ... @@ -1198,40 +1229,85 @@ void GCTAResponseIrf::load_psf(const std::string& filename)
1198 1229 if (gammalib::is_fits(filename)) {
1199 1230  
1200 1231 // Open FITS file
1201   - GFits file(filename);
  1232 + GFits file(fname.filename());
  1233 +
  1234 + // Get the extension name. If an extension name has been specified
  1235 + // then use this name, otherwise use either the
  1236 + // "POINT SPREAD FUNCTION" or the "PSF" extension.
  1237 + std::string extname = "";
  1238 + if (fname.has_extname()) {
  1239 + extname = fname.extname();
  1240 + }
  1241 + else {
  1242 + if (file.contains("POINT SPREAD FUNCTION")) {
  1243 + extname = "POINT SPREAD FUNCTION";
  1244 + }
  1245 + else if (file.contains("PSF")) {
  1246 + extname = "PSF";
  1247 + }
  1248 + }
  1249 +
  1250 + // Continue only if extension name is not empty
  1251 + if (!extname.empty()) {
1202 1252  
1203   - // If file contains a "POINT SPREAD FUNCTION" extension then load it
1204   - // as either a King profile PSF or a 2D PSF
1205   - if (file.contains("POINT SPREAD FUNCTION")) {
1206   - const GFitsTable& table = *file.table("POINT SPREAD FUNCTION");
  1253 + // Get FITS table
  1254 + const GFitsTable& table = *file.table(extname);
  1255 +
  1256 + // Check for King profile specific table columns
1207 1257 if (table.contains("GAMMA") && table.contains("SIGMA")) {
  1258 +
  1259 + // Close FITS file
1208 1260 file.close();
  1261 +
  1262 + // Allocate King profile PSF
1209 1263 m_psf = new GCTAPsfKing(filename);
  1264 +
1210 1265 }
1211   - else if (table.contains("SCALE") && table.contains("SIGMA_1") &&
  1266 +
  1267 + // ... otherwise check for Gaussian profile specific table
  1268 + // columns
  1269 + else if (table.contains("SCALE") && table.contains("SIGMA_1") &&
1212 1270 table.contains("AMPL_2") && table.contains("SIGMA_2") &&
1213 1271 table.contains("AMPL_3") && table.contains("SIGMA_3")) {
  1272 +
  1273 + // Close FITS file
1214 1274 file.close();
  1275 +
  1276 + // Allocate Gaussian profile PSF
1215 1277 m_psf = new GCTAPsf2D(filename);
  1278 +
1216 1279 }
  1280 +
  1281 + // ... otherwise try opening as vector PSF
1217 1282 else {
  1283 +
  1284 + // Close FITS file
1218 1285 file.close();
  1286 +
  1287 + // Allocate vector PSF
  1288 + m_psf = new GCTAPsfVector(filename);
  1289 +
1219 1290 }
1220   - }
1221 1291  
1222   - // ... else load it has PSF vector
  1292 + } // endif: extension name is not empty
  1293 +
  1294 + // Signal that no extension was found
1223 1295 else {
1224   - file.close();
1225   - m_psf = new GCTAPsfVector(filename);
  1296 + std::string msg = "FITS file \""+filename+"\" does not "
  1297 + "contain a valid point spread function table. "
  1298 + "Please specify a valid point spread function "
  1299 + "response file.";
  1300 + throw GException::file_error(G_LOAD_PSF, msg);
1226 1301 }
1227 1302  
1228   - }
  1303 + } // endif: file was FITS file
1229 1304  
1230 1305 // ... otherwise load file as a performance table
1231 1306 else {
1232 1307  
1233 1308 // Allocate a performance table
1234 1309 m_psf = new GCTAPsfPerfTable(filename);
  1310 +
1235 1311 }
1236 1312  
1237 1313 // Return
... ... @@ -1245,16 +1321,19 @@ void GCTAResponseIrf::load_psf(const std::string& filename)
1245 1321 * @param[in] filename Energy dispersion file name.
1246 1322 *
1247 1323 * @exception GException::file_error
1248   - * File not found.
  1324 + * File or extension not found.
1249 1325 *
1250   - * Loads the energy dispersion from a response file. If the file is a FITS
1251   - * file, the method first checks whether it contains a "ENERGY DISPERSION"
1252   - * extension. If such an extension is found, a GCTAEdisp2D energy dispersion
1253   - * is allocated on the basis of the information found in the file. If the
1254   - * "ENERGY DISPERSION" extension is not found, the FITS file will be
1255   - * interpreted as a GCTAEdispRmf energy dispersion. If the file is not a
1256   - * FITS file, it will be interpreted as a GCTAEdispPerfTable performance
1257   - * table.
  1326 + * Loads the energy dispersion from a response file.
  1327 + *
  1328 + * If the file is a FITS file, the method will either use the extension name
  1329 + * specified with the filename, or if no extension name is given, search for
  1330 + * an "ENERGY DISPERSION" or "MATRIX" extension in the file and open the
  1331 + * corresponding FITS table. If columns named "MIGRA_LO" and "MIGRA_HI" are
  1332 + * found in the table, a GCTAEdisp2D object will be allocated. Otherwise, a
  1333 + * GCTAEdispRmf object will be allocated.
  1334 + *
  1335 + * If the file is not a FITS file, it will be interpreted as a
  1336 + * GCTAEdispPerfTable performance table.
1258 1337 ***************************************************************************/
1259 1338 void GCTAResponseIrf::load_edisp(const std::string& filename)
1260 1339 {
... ... @@ -1262,8 +1341,11 @@ void GCTAResponseIrf::load_edisp(const std::string& filename)
1262 1341 if (m_edisp != NULL) delete m_edisp;
1263 1342 m_edisp = NULL;
1264 1343  
  1344 + // Initialise filename
  1345 + GFilename fname(filename);
  1346 +
1265 1347 // Check for existence of file
1266   - if (!gammalib::file_exists_gzip(filename)) {
  1348 + if (!gammalib::file_exists_gzip(fname.filename())) {
1267 1349 std::string msg = "File \""+filename+"\" not found. Please specify "
1268 1350 "a valid energy dispersion response file.";
1269 1351 throw GException::file_error(G_LOAD_EDISP, msg);
... ... @@ -1273,22 +1355,64 @@ void GCTAResponseIrf::load_edisp(const std::string& filename)
1273 1355 if (gammalib::is_fits(filename)) {
1274 1356  
1275 1357 // Open FITS file
1276   - GFits file(filename);
1277   -
1278   - // If file contains an "ENERGY DISPERSION" extension then load it
1279   - // as CTA response table
1280   - if (file.contains("ENERGY DISPERSION")) {
1281   - file.close();
1282   - m_edisp = new GCTAEdisp2D(filename);
  1358 + GFits file(fname.filename());
  1359 +
  1360 + // Get the extension name. If an extension name has been specified
  1361 + // then use this name, otherwise use either the "EFFECTIVE AREA"
  1362 + // or the "SPECRESP" extension.
  1363 + std::string extname = "";
  1364 + if (fname.has_extname()) {
  1365 + extname = fname.extname();
  1366 + }
  1367 + else {
  1368 + if (file.contains("ENERGY DISPERSION")) {
  1369 + extname = "ENERGY DISPERSION";
  1370 + }
  1371 + else if (file.contains("MATRIX")) {
  1372 + extname = "MATRIX";
  1373 + }
1283 1374 }
1284 1375  
1285   - // ... otherwise load it as RMF
  1376 + // Continue only if extension name is not empty
  1377 + if (!extname.empty()) {
  1378 +
  1379 + // Get FITS table
  1380 + const GFitsTable& table = *file.table(extname);
  1381 +
  1382 + // Check for 2D migration matrix
  1383 + if (table.contains("MIGRA_LO") && table.contains("MIGRA_HI")) {
  1384 +
  1385 + // Close FITS file
  1386 + file.close();
  1387 +
  1388 + // Allocate 2D migration matrix
  1389 + m_edisp = new GCTAEdisp2D(filename);
  1390 +
  1391 + }
  1392 +
  1393 + // ... otherwise allocate RMF
  1394 + else {
  1395 +
  1396 + // Close FITS file
  1397 + file.close();
  1398 +
  1399 + // Allocate Gaussian profile PSF
  1400 + m_edisp = new GCTAEdispRmf(filename);
  1401 +
  1402 + }
  1403 +
  1404 + } // endif: extension name is not empty
  1405 +
  1406 + // Signal that no extension was found
1286 1407 else {
1287   - file.close();
1288   - m_edisp = new GCTAEdispRmf(filename);
  1408 + std::string msg = "FITS file \""+filename+"\" does not "
  1409 + "contain a valid energy dispersion table. "
  1410 + "Please specify a valid energy dispersion "
  1411 + "response file.";
  1412 + throw GException::file_error(G_LOAD_EDISP, msg);
1289 1413 }
1290 1414  
1291   - }
  1415 + } // endif: file was FITS file
1292 1416  
1293 1417 // ... otherwise load file as a performance table
1294 1418 else {
... ... @@ -1311,10 +1435,13 @@ void GCTAResponseIrf::load_edisp(const std::string& filename)
1311 1435 * @exception GException::file_error
1312 1436 * File not found.
1313 1437 *
1314   - * Loads the background model from a response file. If the file is a FITS
1315   - * file the method allocates a GCTABackground3D response and loads the
1316   - * information from the FITS file. If the file is not a FITS file it will
1317   - * be interpreted as a GCTABackgroundPerfTable performance table.
  1438 + * Loads the background model from a response file.
  1439 + *
  1440 + * If the file is a FITS file the method allocates a GCTABackground3D
  1441 + * response and loads the information from the FITS file.
  1442 + *
  1443 + * If the file is not a FITS file it will be interpreted as a
  1444 + * GCTABackgroundPerfTable performance table.
1318 1445 ***************************************************************************/
1319 1446 void GCTAResponseIrf::load_background(const std::string& filename)
1320 1447 {
... ... @@ -1322,8 +1449,11 @@ void GCTAResponseIrf::load_background(const std::string& filename)
1322 1449 if (m_background != NULL) delete m_background;
1323 1450 m_background = NULL;
1324 1451  
  1452 + // Initialise filename
  1453 + GFilename fname(filename);
  1454 +
1325 1455 // Check for existence of file
1326   - if (!gammalib::file_exists_gzip(filename)) {
  1456 + if (!gammalib::file_exists_gzip(fname.filename())) {
1327 1457 std::string msg = "File \""+filename+"\" not found. Please specify "
1328 1458 "a valid background response file.";
1329 1459 throw GException::file_error(G_LOAD_BACKGROUND, msg);
... ...