Pierrick Martin / gammalib

Commit fbfd864add43a1c9747cf475c7066e57d985b95f

Authored by Jürgen Knödlseder
1 parent 93799538

Implement "DATE" method in GSPIModelDataSpace (#745)

Inline Side-by-side
Showing 3 changed files with 174 additions and 38 deletions
inst/spi/include/GSPIModelDataSpace.hpp
  View file @fbfd864
... ... @@ -89,45 +89,47 @@ public:
89 89  
90 90 protected:
91 91 // Protected methods
92   - void init_members(void);
93   - void copy_members(const GSPIModelDataSpace& model);
94   - void free_members(void);
95   - void set_pointers(void);
96   - void setup_model(const GObservation& obs) const;
97   - void setup_pars(GSPIEventCube* cube);
98   - void setup_pointing_indices(GSPIEventCube* cube,
  92 + void init_members(void);
  93 + void copy_members(const GSPIModelDataSpace& model);
  94 + void free_members(void);
  95 + void set_pointers(void);
  96 + void setup_model(const GObservation& obs) const;
  97 + void setup_pars(GSPIEventCube* cube);
  98 + void setup_pointing_indices(GSPIEventCube* cube,
  99 + std::vector<int>* indices,
  100 + std::vector<std::string>* names);
  101 + void setup_detector_indices(GSPIEventCube* cube,
  102 + std::vector<int>* indices,
  103 + std::vector<std::string>* names);
  104 + void setup_energy_indices(GSPIEventCube* cube,
99 105 std::vector<int>* indices,
100 106 std::vector<std::string>* names);
101   - void setup_detector_indices(GSPIEventCube* cube,
102   - std::vector<int>* indices,
103   - std::vector<std::string>* names);
104   - void setup_energy_indices(GSPIEventCube* cube,
105   - std::vector<int>* indices,
106   - std::vector<std::string>* names);
107   - void setup_point(GSPIEventCube* cube,
108   - std::vector<int>* indices,
109   - std::vector<std::string>* names);
110   - void setup_orbit(GSPIEventCube* cube,
111   - std::vector<int>* indices,
112   - std::vector<std::string>* names);
113   - void setup_date(GSPIEventCube* cube,
114   - std::vector<int>* indices,
115   - std::vector<std::string>* names);
116   - void add_gedfail(GSPIEventCube* cube,
117   - std::vector<int>* indices,
118   - std::vector<std::string>* names);
119   - void add_gedanneal(GSPIEventCube* cube,
  107 + void setup_point(GSPIEventCube* cube,
  108 + std::vector<int>* indices,
  109 + std::vector<std::string>* names);
  110 + void setup_orbit(GSPIEventCube* cube,
  111 + std::vector<int>* indices,
  112 + std::vector<std::string>* names);
  113 + void setup_date(GSPIEventCube* cube,
  114 + std::vector<int>* indices,
  115 + std::vector<std::string>* names,
  116 + const double& time);
  117 + void add_gedfail(GSPIEventCube* cube,
120 118 std::vector<int>* indices,
121 119 std::vector<std::string>* names);
122   - void setup_dete(GSPIEventCube* cube,
123   - std::vector<int>* indices,
124   - std::vector<std::string>* names);
125   - void setup_evtclass(GSPIEventCube* cube,
126   - std::vector<int>* indices,
127   - std::vector<std::string>* names);
128   - void setup_ebin(GSPIEventCube* cube,
129   - std::vector<int>* indices,
130   - std::vector<std::string>* names);
  120 + void add_gedanneal(GSPIEventCube* cube,
  121 + std::vector<int>* indices,
  122 + std::vector<std::string>* names);
  123 + void setup_dete(GSPIEventCube* cube,
  124 + std::vector<int>* indices,
  125 + std::vector<std::string>* names);
  126 + void setup_evtclass(GSPIEventCube* cube,
  127 + std::vector<int>* indices,
  128 + std::vector<std::string>* names);
  129 + void setup_ebin(GSPIEventCube* cube,
  130 + std::vector<int>* indices,
  131 + std::vector<std::string>* names);
  132 + double get_date_time(const std::string& method) const;
131 133  
132 134 // Protected data members
133 135 mutable GSPIObservation* m_obs; //!< SPI observation
... ...
inst/spi/src/GSPIModelDataSpace.cpp
  View file @fbfd864
... ... @@ -50,6 +50,7 @@ const GModelRegistry g_spi_data_space_registry(&amp;g_spi_data_space_seed);
50 50 #define G_READ "GSPIModelDataSpace::read(GXmlElement&)"
51 51 #define G_WRITE "GSPIModelDataSpace::write(GXmlElement&)"
52 52 #define G_SETUP_MODEL "GSPIModelDataSpace::setup_model(Observation&)"
  53 +#define G_GET_DATE_TIME "GSPIModelDataSpace::get_date_time(std::string&)"
53 54  
54 55 /* __ Macros _____________________________________________________________ */
55 56  
... ... @@ -1017,7 +1018,8 @@ void GSPIModelDataSpace::setup_pointing_indices(GSPIEventCube* cube,
1017 1018  
1018 1019 // Search for "date"
1019 1020 else if (gammalib::contains(method, "date")) {
1020   - setup_date(cube, indices, names);
  1021 + double time = get_date_time(method);
  1022 + setup_date(cube, indices, names, time);
1021 1023 }
1022 1024  
1023 1025 // Handle "gedfail"
... ... @@ -1181,14 +1183,73 @@ void GSPIModelDataSpace::setup_orbit(GSPIEventCube* cube,
1181 1183 * @param[in] cube Event cube.
1182 1184 * @param[in,out] indices Vector of pointing indices.
1183 1185 * @param[in,out] names Vector of pointing names.
  1186 + * @param[in] time Time scale in seconds.
1184 1187 *
1185 1188 * Setup vectors of pointing indices and names for "date" method.
1186 1189 ***************************************************************************/
1187 1190 void GSPIModelDataSpace::setup_date(GSPIEventCube* cube,
1188 1191 std::vector<int>* indices,
1189   - std::vector<std::string>* names)
  1192 + std::vector<std::string>* names,
  1193 + const double& time)
1190 1194 {
1191   - // TODO: implement "date" method
  1195 + // Compute start time, number of time bins and length of time bin in
  1196 + // seconds
  1197 + double tstart = cube->gti().tstart().secs();
  1198 + double tbin = cube->gti().telapse();
  1199 + int nbins = int(tbin / time + 0.5);
  1200 + if (nbins < 1) {
  1201 + nbins = 1;
  1202 + }
  1203 + tbin /= double(nbins);
  1204 +
  1205 + // Get number of pointings
  1206 + int npt = indices->size();
  1207 +
  1208 + // Loop over all pointings
  1209 + for (int ipt = 0; ipt < npt; ++ipt) {
  1210 +
  1211 + // Get time of pointing in seconds
  1212 + double t = 0.5 * (cube->gti().tstart(ipt).secs() +
  1213 + cube->gti().tstop(ipt).secs());
  1214 +
  1215 + // Compute grouping index. Make sure that it is in valid range.
  1216 + int index = int((t - tstart) / tbin);
  1217 + if (index < 0) {
  1218 + index = 0;
  1219 + }
  1220 + else if (index >= nbins) {
  1221 + index = nbins - 1;
  1222 + }
  1223 +
  1224 + // Store index
  1225 + (*indices)[ipt] = index;
  1226 +
  1227 + } // endfor: looped over indices
  1228 +
  1229 + // Remove unused indices
  1230 + for (int index = 0, used_index = 0; index < nbins; ++index) {
  1231 +
  1232 + // Replace all indices "index" by "used_index"
  1233 + int nindex = 0;
  1234 + for (int ipt = 0; ipt < npt; ++ipt) {
  1235 + if ((*indices)[ipt] == index) {
  1236 + (*indices)[ipt] = used_index;
  1237 + nindex++;
  1238 + }
  1239 + }
  1240 +
  1241 + // If we found indices to replaced the increase "used_index"
  1242 + if (nindex > 0) {
  1243 +
  1244 + // Add date to names
  1245 + names->push_back("T" + gammalib::str(used_index, "%6.6d"));
  1246 +
  1247 + // Increment used index
  1248 + used_index++;
  1249 +
  1250 + } // endif: index was used
  1251 +
  1252 + } // endfor: looped over all pointings
1192 1253  
1193 1254 // Return
1194 1255 return;
... ... @@ -1309,3 +1370,56 @@ void GSPIModelDataSpace::setup_ebin(GSPIEventCube* cube,
1309 1370 // Return
1310 1371 return;
1311 1372 }
  1373 +
  1374 +
  1375 +/***********************************************************************//**
  1376 + * @brief Get time scale from method string
  1377 + *
  1378 + * @param[in] method Method string.
  1379 + * @return Time scale (seconds).
  1380 + *
  1381 + * Get the time scale for the "DATE" method from the method string. The
  1382 + * following time units are supported: min, hour, day, week, month or year.
  1383 + ***************************************************************************/
  1384 +double GSPIModelDataSpace::get_date_time(const std::string& method) const
  1385 +{
  1386 + // Get position of "date" method in string
  1387 + size_t start = method.find("date") + 4;
  1388 +
  1389 + // Get scale of time
  1390 + double scale = 0.0;
  1391 + size_t stop = std::string::npos;
  1392 + if ((stop = method.find("min", start)) != std::string::npos) {
  1393 + scale = 60.0;
  1394 + }
  1395 + else if ((stop = method.find("hour", start)) != std::string::npos) {
  1396 + scale = 3600.0;
  1397 + }
  1398 + else if ((stop = method.find("day", start)) != std::string::npos) {
  1399 + scale = 86400.0;
  1400 + }
  1401 + else if ((stop = method.find("week", start)) != std::string::npos) {
  1402 + scale = 604800.0;
  1403 + }
  1404 + else if ((stop = method.find("month", start)) != std::string::npos) {
  1405 + scale = 2628000.0;
  1406 + }
  1407 + else if ((stop = method.find("year", start)) != std::string::npos) {
  1408 + scale = 31536000.0;
  1409 + }
  1410 +
  1411 + // Throw an exception if no time scale was found
  1412 + if (stop == std::string::npos) {
  1413 + std::string msg = "Method string \""+method+"\" does not contain "
  1414 + "time units for \"DATE\" method. Please specify one "
  1415 + "of the following units: min, hour, day, week, "
  1416 + "month or year.";
  1417 + throw GException::invalid_value(G_GET_DATE_TIME, msg);
  1418 + }
  1419 +
  1420 + // Get time in seconds
  1421 + double time = gammalib::todouble(method.substr(start, stop-start)) * scale;
  1422 +
  1423 + // Return time
  1424 + return time;
  1425 +}
... ...
inst/spi/test/test_SPI.cpp
  View file @fbfd864
... ... @@ -271,6 +271,26 @@ void TestGSPI::test_modeldataspace(void)
271 271 test_value(model6.size(), 88, "POINT instance size");
272 272 test_value(model6[35].name(), "GEDSAT P000035", "POINT instance parameter");
273 273  
  274 + // Construct "date 10 min" model from OG
  275 + GSPIModelDataSpace model7(obs, "GEDSAT", "date 10 mins", 0);
  276 + test_value(model7.size(), 88, "DATE 10 MIN instance size");
  277 + test_value(model7[35].name(), "GEDSAT T000035", "DATE 10 MIN instance parameter");
  278 +
  279 + // Construct "date 3 hours" model from OG
  280 + GSPIModelDataSpace model8(obs, "GEDSAT", "date 3 hours", 0);
  281 + test_value(model8.size(), 20, "DATE 3 HOURS instance size");
  282 + test_value(model8[5].name(), "GEDSAT T000005", "DATE 3 HOURS instance parameter");
  283 +
  284 + // Construct "date 1 day" model from OG
  285 + GSPIModelDataSpace model9(obs, "GEDSAT", "date 1 day", 0);
  286 + test_value(model9.size(), 3, "DATE 1 DAY instance size");
  287 + test_value(model9[1].name(), "GEDSAT T000001", "DATE 1 DAY instance parameter");
  288 +
  289 + // Construct "date 2 weeks" model from OG
  290 + GSPIModelDataSpace model10(obs, "GEDSAT", "date 2 weeks", 0);
  291 + test_value(model10.size(), 1, "DATE 2 WEEKS instance size");
  292 + test_value(model10[0].name(), "GEDSAT T000000", "DATE 2 WEEKS instance parameter");
  293 +
274 294 // Return
275 295 return;
276 296 }
... ...