Justine Devin / gammalib

Commit e6e0f0dc6add04bf6900ebb700154c84f4d55ddd

Authored by knodlseder
1 parent aaf1d231

add read_ds_roi() and read_ds_ebounds() methods and read DS keywords in unbinned loading

Inline Side-by-side
Showing 1 changed file with 188 additions and 3 deletions
inst/cta/src/GCTAObservation.cpp
  View file @e6e0f0d
@@ -33,6 +33,8 @@ @@ -33,6 +33,8 @@
33 #include "GIntegrand.hpp" 33 #include "GIntegrand.hpp"
34 34
35 /* __ Method name definitions ____________________________________________ */ 35 /* __ Method name definitions ____________________________________________ */
  36 +#define G_READ_DS_EBOUNDS "GCTAObservation::read_ds_ebounds(GFitsHDU*)"
  37 +#define G_READ_DS_ROI "GCTAObservation::read_ds_roi(GFitsHDU*)"
36 38
37 /* __ Macros _____________________________________________________________ */ 39 /* __ Macros _____________________________________________________________ */
38 40
@@ -310,9 +312,15 @@ std::string GCTAObservation::print(void) const @@ -310,9 +312,15 @@ std::string GCTAObservation::print(void) const
310 ***************************************************************************/ 312 ***************************************************************************/
311 void GCTAObservation::load_unbinned(const std::string& filename) 313 void GCTAObservation::load_unbinned(const std::string& filename)
312 { 314 {
313 - // Delete old events. We do not call clear() here since we want to  
314 - // preserve any existing response function. 315 + // Delete old events and ROI. We do not call clear() here since we want
  316 + // to preserve any existing response function.
315 if (m_events != NULL) delete m_events; 317 if (m_events != NULL) delete m_events;
  318 + if (m_roi != NULL) delete m_roi;
  319 + m_events = NULL;
  320 + m_roi = NULL;
  321 +
  322 + // Reset energy boundaries
  323 + m_ebounds.clear();
316 324
317 // Allocate events 325 // Allocate events
318 GCTAEventList* events = new GCTAEventList; 326 GCTAEventList* events = new GCTAEventList;
@@ -329,7 +337,11 @@ void GCTAObservation::load_unbinned(const std::string& filename) @@ -329,7 +337,11 @@ void GCTAObservation::load_unbinned(const std::string& filename)
329 337
330 // Read observation attributes 338 // Read observation attributes
331 read_attributes(hdu); 339 read_attributes(hdu);
332 - 340 +
  341 + // Read data selection keywords
  342 + read_ds_ebounds(hdu);
  343 + read_ds_roi(hdu);
  344 +
333 // Close FITS file 345 // Close FITS file
334 file.close(); 346 file.close();
335 347
@@ -485,6 +497,8 @@ void GCTAObservation::free_members(void) @@ -485,6 +497,8 @@ void GCTAObservation::free_members(void)
485 * @brief Read observation attributes 497 * @brief Read observation attributes
486 * 498 *
487 * @param[in] hdu FITS HDU 499 * @param[in] hdu FITS HDU
  500 + *
  501 + * Reads the observation attributes.
488 ***************************************************************************/ 502 ***************************************************************************/
489 void GCTAObservation::read_attributes(GFitsHDU* hdu) 503 void GCTAObservation::read_attributes(GFitsHDU* hdu)
490 { 504 {
@@ -519,6 +533,177 @@ void GCTAObservation::read_attributes(GFitsHDU* hdu) @@ -519,6 +533,177 @@ void GCTAObservation::read_attributes(GFitsHDU* hdu)
519 533
520 534
521 /***********************************************************************//** 535 /***********************************************************************//**
  536 + * @brief Read energy boundary data selection keywords
  537 + *
  538 + * @param[in] hdu FITS HDU
  539 + *
  540 + * @exception GCTAException::no_ebds
  541 + * Invalid energy data selection encountered
  542 + *
  543 + * Reads the energy boundary data selection keywords by searching for a
  544 + * DSTYPx keyword named "ENERGY". The data selection information is expected
  545 + * to be in the format "200:50000", where the 2 arguments are the minimum
  546 + * and maximum energy. The energy unit is given by the keyword DSUNIx, which
  547 + * supports keV, MeV, GeV and TeV (case independent). No detailed syntax
  548 + * checking is performed.
  549 + *
  550 + * This method does nothing if the HDU is NULL.
  551 + ***************************************************************************/
  552 +void GCTAObservation::read_ds_ebounds(GFitsHDU* hdu)
  553 +{
  554 + // Continue only if HDU is valid
  555 + if (hdu != NULL) {
  556 +
  557 + // Reset energy boundaries
  558 + m_ebounds.clear();
  559 +
  560 + // Get number of data selection keywords
  561 + int ndskeys = 0;
  562 + try {
  563 + ndskeys = hdu->integer("NDSKEYS");
  564 + }
  565 + catch (GException::fits_key_not_found) {
  566 + ;
  567 + }
  568 +
  569 + // Loop over all data selection keys
  570 + for (int i = 1; i <= ndskeys; ++i) {
  571 + std::string type_key = "DSTYP"+str(i);
  572 + std::string unit_key = "DSUNI"+str(i);
  573 + std::string value_key = "DSVAL"+str(i);
  574 + try {
  575 + if (hdu->string(type_key) == "ENERGY") {
  576 + std::string unit = toupper(hdu->string(unit_key));
  577 + std::string value = hdu->string(value_key);
  578 + std::vector<std::string> ebounds = split(value, ":");
  579 + if (ebounds.size() == 2) {
  580 + double emin = todouble(ebounds[0]);
  581 + double emax = todouble(ebounds[1]);
  582 + GEnergy e_min;
  583 + GEnergy e_max;
  584 + if (unit == "KEV") {
  585 + e_min.keV(emin);
  586 + e_max.keV(emax);
  587 + }
  588 + else if (unit == "MEV") {
  589 + e_min.MeV(emin);
  590 + e_max.MeV(emax);
  591 + }
  592 + else if (unit == "GEV") {
  593 + e_min.GeV(emin);
  594 + e_max.GeV(emax);
  595 + }
  596 + else if (unit == "TEV") {
  597 + e_min.TeV(emin);
  598 + e_max.TeV(emax);
  599 + }
  600 + else {
  601 + throw GCTAException::no_ebds(G_READ_DS_EBOUNDS,
  602 + "Invalid energy unit \""+unit+
  603 + "\" encountered in data selection key \""+
  604 + unit_key+"\"");
  605 + }
  606 + m_ebounds.append(e_min, e_max);
  607 + }
  608 + else {
  609 + throw GCTAException::no_ebds(G_READ_DS_EBOUNDS,
  610 + "Invalid energy value \""+value+
  611 + "\" encountered in data selection key \""+
  612 + value_key+"\"");
  613 + }
  614 + } // endif: ENERGY type found
  615 + }
  616 + catch (GException::fits_key_not_found) {
  617 + ;
  618 + }
  619 + } // endfor: looped over data selection keys
  620 +
  621 + } // endif: HDU was valid
  622 +
  623 + // Return
  624 + return;
  625 +}
  626 +
  627 +
  628 +/***********************************************************************//**
  629 + * @brief Read ROI data selection keywords
  630 + *
  631 + * @param[in] hdu FITS HDU
  632 + *
  633 + * @exception GException::no_roi
  634 + * Invalid ROI data selection encountered
  635 + *
  636 + * Reads the ROI data selection keywords by searching for a DSTYPx keyword
  637 + * named "POS(RA,DEC)". The data selection information is expected to be
  638 + * in the format "CIRCLE(267.0208,-24.78,4.5)", where the 3 arguments are
  639 + * Right Ascension, Declination and radius in units of degrees. No detailed
  640 + * syntax checking is performed.
  641 + *
  642 + * This method does nothing if the HDU is NULL.
  643 + ***************************************************************************/
  644 +void GCTAObservation::read_ds_roi(GFitsHDU* hdu)
  645 +{
  646 + // Continue only if HDU is valid
  647 + if (hdu != NULL) {
  648 +
  649 + // Delete any old ROI
  650 + if (m_roi != NULL) delete m_roi;
  651 +
  652 + // Allocate ROI
  653 + GCTARoi* roi = new GCTARoi;
  654 + m_roi = roi;
  655 +
  656 + // Get number of data selection keywords
  657 + int ndskeys = 0;
  658 + try {
  659 + ndskeys = hdu->integer("NDSKEYS");
  660 + }
  661 + catch (GException::fits_key_not_found) {
  662 + ;
  663 + }
  664 +
  665 + // Loop over all data selection keys
  666 + for (int i = 1; i <= ndskeys; ++i) {
  667 + std::string type_key = "DSTYP"+str(i);
  668 + std::string unit_key = "DSUNI"+str(i);
  669 + std::string value_key = "DSVAL"+str(i);
  670 + try {
  671 + if (hdu->string(type_key) == "POS(RA,DEC)") {
  672 + std::string unit = toupper(hdu->string(unit_key));
  673 + std::string value = hdu->string(value_key);
  674 + value = strip_chars(value, "CIRCLE(");
  675 + value = strip_chars(value, ")");
  676 + std::vector<std::string> args = split(value, ",");
  677 + if (args.size() == 3) {
  678 + double ra = todouble(args[0]);
  679 + double dec = todouble(args[1]);
  680 + double rad = todouble(args[2]);
  681 + GCTAInstDir dir;
  682 + dir.radec_deg(ra, dec);
  683 + roi->centre(dir);
  684 + roi->radius(rad);
  685 + }
  686 + else {
  687 + throw GException::no_roi(G_READ_DS_ROI,
  688 + "Invalid acceptance cone value \""+value+
  689 + "\" encountered in data selection key \""+
  690 + value_key+"\"");
  691 + }
  692 + } // endif: POS(RA,DEC) type found
  693 + }
  694 + catch (GException::fits_key_not_found) {
  695 + ;
  696 + }
  697 + } // endfor: looped over data selection keys
  698 +
  699 + } // endif: HDU was valid
  700 +
  701 + // Return
  702 + return;
  703 +}
  704 +
  705 +
  706 +/***********************************************************************//**
522 * @brief Write observation attributes 707 * @brief Write observation attributes
523 * 708 *
524 * @param[in] hdu FITS HDU 709 * @param[in] hdu FITS HDU