GCOSEventList.cpp

Go to the documentation of this file.

/***************************************************************************
 *                GCOSEventList.cpp - COSI event list class                *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 2026 by Juergen Knoedlseder                              *
 * ----------------------------------------------------------------------- *
 *                                                                         *
 *  This program is free software: you can redistribute it and/or modify   *
 *  it under the terms of the GNU General Public License as published by   *
 *  the Free Software Foundation, either version 3 of the License, or      *
 *  (at your option) any later version.                                    *
 *                                                                         *
 *  This program is distributed in the hope that it will be useful,        *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 *  GNU General Public License for more details.                           *
 *                                                                         *
 *  You should have received a copy of the GNU General Public License      *
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
 *                                                                         *
 ***************************************************************************/
/**
 * @file GCOSEventList.hpp
 * @brief COSI event list class implementation
 * @author Juergen Knoedlseder
 */
 
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <algorithm>   // std::sort
#include <typeinfo>
#include "GFits.hpp"
#include "GException.hpp"
#include "GFitsBinTable.hpp"
#include "GFitsTableDoubleCol.hpp"
#include "GCOSTools.hpp"
#include "GCOSEventList.hpp"
 
/* __ Method name definitions ____________________________________________ */
#define G_OPERATOR                          "GCOSEventList::operator[](int&)"
#define G_READ                                  "GCOSEventList::read(GFits&)"
#define G_ROI                                     "GCOSEventList::roi(GRoi&)"
 
/* __ Macros _____________________________________________________________ */
 
/* __ Coding definitions _________________________________________________ */
 
/* __ Debug definitions __________________________________________________ */
 
 
/*==========================================================================
 =                                                                         =
 =                         Constructors/destructors                        =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Void constructor
 *
 * Creates an empty COSI event list.
 ***************************************************************************/
GCOSEventList::GCOSEventList(void) : GEventList()
{
    // Initialise class members for clean destruction
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief File name constructor
 *
 * @param[in] filename COSI event list filename.
 *
 * Construct COSI event list object by loading the events from a
 * FITS file.
 ***************************************************************************/
GCOSEventList::GCOSEventList(const GFilename& filename) : GEventList()
{
    // Initialise members
    init_members();
 
    // Load event list
    load(filename);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy constructor
 *
 * @param[in] list COSI event list.
 ***************************************************************************/
GCOSEventList::GCOSEventList(const GCOSEventList& list) : GEventList(list)
{
    // Initialise members
    init_members();
 
    // Copy members
    copy_members(list);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Destructor
 ***************************************************************************/
GCOSEventList::~GCOSEventList(void)
{
    // Free members
    free_members();
 
    // Return
    return;
}
 
 
/*==========================================================================
 =                                                                         =
 =                               Operators                                 =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Assignment operator
 *
 * @param[in] list COSI event list.
 * @return COSI event list.
 ***************************************************************************/
GCOSEventList& GCOSEventList::operator=(const GCOSEventList& list)
{
    // Execute only if object is not identical
    if (this != &list) {
 
        // Copy base class members
        this->GEventList::operator=(list);
 
        // Free members
        free_members();
 
        // Initialise members
        init_members();
 
        // Copy members
        copy_members(list);
 
    } // endif: object was not identical
 
    // Return this object
    return *this;
}
 
 
/***********************************************************************//**
 * @brief COSI event atom access operator
 *
 * @param[in] index Event index [0,...,size()-1].
 * @return Pointer to COSI event atom.
 *
 * @exception GException::out_of_range
 *            Event index outside valid range.
 *
 * Returns pointer to a COSI event atom.
 ***************************************************************************/
GCOSEventAtom* GCOSEventList::operator[](const int& index)
{
    // Optionally check if the index is valid
    #if defined(G_RANGE_CHECK)
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_OPERATOR, "Event index", index, size());
    }
    #endif
 
    // Return pointer
    return (&(m_events[index]));
}
 
 
/***********************************************************************//**
 * @brief COSI event atom access operator
 *
 * @param[in] index Event index [0,...,size()-1].
 * @return Pointer to COSI event atom.
 *
 * @exception GException::out_of_range
 *            Event index outside valid range.
 *
 * Returns pointer to a COSI event atom.
 ***************************************************************************/
const GCOSEventAtom* GCOSEventList::operator[](const int& index) const
{
    // Optionally check if the index is valid
    #if defined(G_RANGE_CHECK)
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_OPERATOR, "Event index", index, size());
    }
    #endif
 
    // Return pointer
    return (&(m_events[index]));
}
 
 
/*==========================================================================
 =                                                                         =
 =                             Public methods                              =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Clear COSI event list
 *
 * Clears COSI event list by resetting all class members to an
 * initial state. Any information that was present before will be lost.
 ***************************************************************************/
void GCOSEventList::clear(void)
{
    // Free class members (base and derived classes, derived class first)
    free_members();
    this->GEventList::free_members();
    this->GEvents::free_members();
 
    // Initialise members
    this->GEvents::init_members();
    this->GEventList::init_members();
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Clone event list
 *
 * @return Pointer to deep copy of COSI event list.
 ***************************************************************************/
GCOSEventList* GCOSEventList::clone(void) const
{
    return new GCOSEventList(*this);
}
 
 
/***********************************************************************//**
 * @brief Load COSI events from FITS file
 *
 * @param[in] filename COSI event list FITS file name.
 *
 * Loads COSI events from a FITS file into the event list.
 ***************************************************************************/
void GCOSEventList::load(const GFilename& filename)
{
    // Open FITS file
    GFits fits(filename);
 
    // Read event list from FITS file
    read(fits);
 
    // Close FITS file
    fits.close();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Save COSI events
 *
 * @param[in] filename COSI event list FITS file name.
 * @param[in] clobber Overwrite existing FITS file?
 *
 * Save COSI events into a FITS file.
 ***************************************************************************/
void GCOSEventList::save(const GFilename& filename,
                         const bool&      clobber) const
{
    // Open FITS file
    GFits fits;
 
    // Write events
    write(fits);
 
    // Save events
    fits.saveto(filename, clobber);
 
    // Close FITS file
    fits.close();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read COSI events from FITS file.
 *
 * @param[in] fits FITS file.
 *
 * @exception GException::invalid_value
 *            Invalid data space keyword formatting encountered
 *
 * Read the COSI event list from a FITS file object.
 ***************************************************************************/
void GCOSEventList::read(const GFits& fits)
{
    // Clear object
    clear();
 
    // Initialise energy boundaries and RoI radius
    GEnergy emin;
    GEnergy emax;
    double  radius = 180.0;
 
    // Get reference to event table
    const GFitsTable& table = *fits.table(1);
 
    // Extract number of events in FITS file
    int num = table.nrows();
 
    // If there are events then read them
    if (num > 0) {
 
        // Initialise boundaries
        double phi_min = 0.0;
        double phi_max = 0.0;
 
        // Reserve data
        m_events.reserve(num);
 
        // Get column pointers
        const GFitsTableCol* ptr_energies  = table["Energies"];     // keV
        const GFitsTableCol* ptr_time      = table["TimeTags"];     // sec
        const GFitsTableCol* ptr_phi       = table["Phi"];          // rad
        const GFitsTableCol* ptr_chi_local = table["Chi local"];    // rad
        const GFitsTableCol* ptr_psi_local = table["Psi local"];    // rad
        const GFitsTableCol* ptr_chi_gal   = table["Chi galactic"]; // deg
        const GFitsTableCol* ptr_psi_gal   = table["Psi galactic"]; // deg
 
        // Copy data from columns into GCOSEventAtom objects
        for (int i = 0; i < num; ++i) {
 
            // Set instrument direction
            GSkyDir dir_gal;   // Chi,Psi in galactic coordinates
            GSkyDir dir_local; // Chi,Psi in spacecraft coordinates
            double  phi = ptr_phi->real(i) * gammalib::rad2deg;
            dir_gal.lb_deg(ptr_chi_gal->real(i), ptr_psi_gal->real(i));
            dir_local.radec(ptr_chi_local->real(i), ptr_psi_local->real(i));
            GCOSInstDir dir(dir_gal, dir_local, phi);
 
            // Set energy
            GEnergy energy;
            energy.keV(ptr_energies->real(i));
 
            // Set time
            GTime time(gammalib::cos_time(ptr_time->real(i)));
 
            // Allocate event
            GCOSEventAtom event(dir, energy, time);
 
            // Append event
            m_events.push_back(event);
 
            // Update boundaries
            if (i == 0) {
                emin    = energy;
                emax    = energy;
                phi_min = phi;
                phi_max = phi;
            }
            else {
                if (energy < emin) {
                    emin = energy;
                }
                if (energy > emax) {
                    emax = energy;
                }
                if (phi < phi_min) {
                    phi_min = phi;
                }
                if (phi > phi_max) {
                    phi_max = phi;
                }
            }
 
        } // endfor: looped over all events
 
    } // endif: there were events
 
    // If energy selection keywords are present then read energy boundaries
    // from keywords
    if (table.has_card("DSUNI2") && table.has_card("DSVAL2")) {
        std::string unit                = table.string("DSUNI2");
        std::string dsval               = table.string("DSVAL2");
        std::vector<std::string> energy = gammalib::split(dsval, ":");
        if (energy.size() != 2) {
            std::string msg = "Invalid data selecetion value string DSVAL2="+
                              dsval+". Please specify a FITS file with a "
                              "correct formatting.";
            throw GException::invalid_value(G_READ, msg);
        }
        emin = GEnergy(gammalib::todouble(energy[0]), unit);
        emax = GEnergy(gammalib::todouble(energy[1]), unit);
    }
 
    // Set energy boundaries
    m_ebounds = GEbounds(emin, emax);
 
    // If GTI extension is present in FITS file then read Good Time Intervals
    // from that extension
    if (fits.contains(gammalib::extname_gti)) {
        const GFitsTable& gti = *fits.table(gammalib::extname_gti);
        m_gti.read(gti);
    }
 
    // If RoI keywords are present then read radius from keywords
    if (table.has_card("DSVAL3")) {
        std::string dsval               = table.string("DSVAL3");
        dsval                           = gammalib::replace_segment(dsval, "CIRCLE(", "");
        dsval                           = gammalib::rstrip_chars(dsval, ")");
        std::vector<std::string> angles = gammalib::split(dsval, ",");
        if (angles.size() != 3) {
            std::string msg = "Invalid data selecetion value string DSVAL3="+
                              dsval+". Please specify a FITS file with a "
                              "correct formatting.";
            throw GException::invalid_value(G_READ, msg);
        }
        radius = gammalib::todouble(angles[2]);
    }
 
    // Set RoI radius
    m_roi.radius(radius);
 
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Write COSI event list into FITS file.
 *
 * @param[in] fits FITS file.
 *
 * Write the COSI event list into FITS file.
 ***************************************************************************/
void GCOSEventList::write(GFits& fits) const
{
    // Allocate empty FITS binary table
    GFitsBinTable table;
 
    // Set FITS extension name
    table.extname("EVENT LIST");
 
    // Allocate columns
    GFitsTableDoubleCol col_energies("Energies", size());
    GFitsTableDoubleCol col_time("TimeTags", size());
    GFitsTableDoubleCol col_phi("Phi", size());
    GFitsTableDoubleCol col_chi_local("Chi local", size());
    GFitsTableDoubleCol col_psi_local("Psi local", size());
    GFitsTableDoubleCol col_chi_gal("Chi galactic", size());
    GFitsTableDoubleCol col_psi_gal("Psi galactic", size());
 
    // Set units of columns
    // (see http://fits.gsfc.nasa.gov/standard30/fits_standard30aa.pdf)
    col_energies.unit("keV");
    col_time.unit("s");
    col_phi.unit("rad");
    col_chi_local.unit("rad");
    col_psi_local.unit("rad");
    col_chi_gal.unit("deg");
    col_psi_gal.unit("deg");
 
    // Fill columns
    for (int i = 0; i < size(); ++i) {
        col_energies(i)  = m_events[i].energy().keV();
        col_time(i)      = gammalib::cos_seconds(m_events[i].time());
        col_phi(i)       = m_events[i].dir().phi() * gammalib::deg2rad;
        col_chi_local(i) = m_events[i].dir().dir_local().ra();
        col_psi_local(i) = m_events[i].dir().dir_local().dec();
        col_chi_gal(i)   = m_events[i].dir().dir().l_deg();
        col_psi_gal(i)   = m_events[i].dir().dir().b_deg();
    }
 
    // Append columns to table
    table.append(col_energies);
    table.append(col_time);
    table.append(col_phi);
    table.append(col_chi_local);
    table.append(col_psi_local);
    table.append(col_chi_gal);
    table.append(col_psi_gal);
 
    // Set data selection strings
    std::string dsref1 = ":"+gammalib::extname_gti;
    std::string dsval2 = gammalib::str(emin().keV()) + ":" + gammalib::str(emax().keV());
    std::string dsval3 = "CIRCLE(0,0," + gammalib::str(m_roi.radius()) + ")";
 
    // Write data selection strings
    table.card("DSTYP1", "TIME",   "Data sub-space type");
    table.card("DSUNI1", "s",      "Data sub-space unit");
    table.card("DSVAL1", "TABLE",  "Data sub-space value");
    table.card("DSREF1", dsref1,   "Data sub-space reference");
    table.card("DSTYP2", "ENERGY", "Data sub-space type");
    table.card("DSUNI2", "keV",    "Data sub-space unit");
    table.card("DSVAL2", dsval2,   "Data sub-space value");
    table.card("DSTYP3", "POS(RA,DEC)", "Data sub-space type");
    table.card("DSUNI3", "deg",         "Data sub-space unit");
    table.card("DSVAL3", dsval3,        "Data sub-space value");
    table.card("NDSKEYS", 3, "Number of data sub-space keys");
 
    // Append event table to FITS file
    fits.append(table);
 
    // Write GTI extension
    gti().write(fits, gammalib::extname_gti);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Set region of interest
 *
 * @param[in] roi Region of interest.
 *
 * @exception GException::invalid_argument
 *            Specified RoI is not a COSI RoI.
 *
 * Sets the region of interest for the observation.
 ***************************************************************************/
void GCOSEventList::roi(const GRoi& roi)
{
    // Get pointer on COSI region of interest
    const GCOSRoi* cosroi = dynamic_cast<const GCOSRoi*>(&roi);
    if (cosroi == NULL) {
        std::string cls = std::string(typeid(&roi).name());
        std::string msg = "Region of interest of type \""+cls+"\" is "
                          "not a COSI RoI. Please specify a "
                          "COSI RoI as argument.";
        throw GException::invalid_argument(G_ROI, msg);
    }
 
    // Set RoI
    m_roi = *cosroi;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Append event to event list
 *
 * @param[in] event Event.
 *
 * Appends an event to the end of the event list.
 ***************************************************************************/
void GCOSEventList::append(const GCOSEventAtom& event)
{
    // Append event
    m_events.push_back(event);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Extend event list with another event list
 *
 * @param[in] list Event list.
 *
 * Extend event list with another event list. All events of @p list are
 * appended after the existing events. If a sorted event list is required
 * then call the sort() method after this method.
 ***************************************************************************/
void GCOSEventList::extend(const GCOSEventList& list)
{
    // Get number of events in specified event list
    int num = list.size();
 
    // Continue only if there are elements in the event list
    if (num > 0) {
 
        // Reserve enough space
        reserve(size() + num);
 
        // Append events
        for (int i = 0; i < num; ++i) {
            m_events.push_back(*(list[i]));
        }
 
    } // endif: event list was empty
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Remove events from event list
 *
 * @param[in] index Index from which on events should be removed.
 * @param[in] number Number of event to remove (default: 1).
 *
 * Removes events from the event list. This method does nothing if @p index
 * points beyond the event list. The method does also gently accept
 * @p number arguments where @p index + @p number reach beyond the event
 * list. In that case, all events from event @p index on will be removed.
 ***************************************************************************/
void GCOSEventList::remove(const int& index, const int& number)
{
    // Continue only if index is valid
    if (index < size()) {
 
        // Determine number of elements to remove
        int n_remove = (index + number > size()) ? size() - index : number;
 
        // Remove events
        m_events.erase(m_events.begin() + index,
                       m_events.begin() + index + n_remove);
 
    } // endif: index was valid
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Sort events by increasing time
 *
 * Sorts the events in the event list by increasing time value.
 ***************************************************************************/
void GCOSEventList::sort(void)
{
    // Sort events in event list
    std::sort(m_events.begin(), m_events.end());
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Print COSI event list information
 *
 * @param[in] chatter Chattiness.
 * @return String containing COSI event list information.
 ***************************************************************************/
std::string GCOSEventList::print(const GChatter& chatter) const
{
    // Initialise result string
    std::string result;
 
    // Continue only if chatter is not silent
    if (chatter != SILENT) {
 
        // Append header
        result.append("=== GCOSEventList ===");
 
        // Append information
        result.append("\n"+gammalib::parformat("Number of events") +
                      gammalib::str(number()));
 
        // Append time information
        if (gti().size() > 0) {
            result.append("\n"+gammalib::parformat("MJD interval"));
            result.append(gammalib::str(tstart().mjd())+" - ");
            result.append(gammalib::str(tstop().mjd())+" days");
            result.append("\n"+gammalib::parformat("UTC interval"));
            result.append(tstart().utc()+" - ");
            result.append(tstop().utc());
        }
        else {
            result.append("\n"+gammalib::parformat("MJD interval"));
            result.append("not defined");
        }
 
        // Append energy information
        result.append("\n"+gammalib::parformat("Energy interval"));
        if (ebounds().size() > 0) {
            result.append(gammalib::str(ebounds().emin().MeV()));
            result.append(" - ");
            result.append(gammalib::str(ebounds().emax().MeV())+" MeV");
        }
        else {
            result.append("not defined");
        }
 
    } // endif: chatter was not silent
 
    // Return result
    return result;
}
 
 
/*==========================================================================
 =                                                                         =
 =                            Private methods                              =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Initialise class members
 ***************************************************************************/
void GCOSEventList::init_members(void)
{
    // Initialise members
    m_roi.clear();
    m_events.clear();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy class members
 *
 * @param[in] list COSI event list.
 ***************************************************************************/
void GCOSEventList::copy_members(const GCOSEventList& list)
{
    // Copy members
    m_roi    = list.m_roi;
    m_events = list.m_events;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Delete class members
 ***************************************************************************/
void GCOSEventList::free_members(void)
{
    // Return
    return;
}