From 8f3f4932f78b27f75a684c5d8bc3eb8aa37ca6aa Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 14 Apr 2017 14:22:28 +0200 Subject: [PATCH] Add support for Monte Carlo identifier metadata (#1483) --- ChangeLog | 3 ++- NEWS | 7 ++++++- inst/cta/include/GCTAEventList.hpp | 18 ++++++++++++------ inst/cta/pyext/GCTAEventList.i | 2 ++ inst/cta/src/GCTAEventAtom.cpp | 4 +++- inst/cta/src/GCTAEventList.cpp | 167 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------- 6 files changed, 181 insertions(+), 20 deletions(-) diff --git a/ChangeLog b/ChangeLog index 5beabda..203e17a 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2017-04-13 +2017-04-14 * Version 1.3.0 released ========================= + Add support for Monte Carlo identifier metadata (#1483) Add GEventList::has_mc_id() and GEventAtom::mc_id() methods (#1483) Add GEventList::has_phase() and GEventList::has_detxy() methods (#1982) Add GEventList::append_column() method (#1991) diff --git a/NEWS b/NEWS index 426c089..9a32f9a 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.3.0.dev1 -13 April 2017 +14 April 2017 1. Introduction @@ -54,6 +54,7 @@ The following methods have been added: - GCTAEventList::has_detxy() - GCTAEventList::has_mc_id(bool&) - GCTAEventList::has_mc_id() +- GCTAEventList::set_mc_id_names(std::vector&, std::vector&) - GCTAEventAtom::mc_id(int&) - GCTAEventAtom::mc_id() @@ -233,6 +234,10 @@ Add GCTAEventList::has_mc_id() and GCTAEventAtom::mc_id() getter and setter meth to set the Monte Carlo identifier of a given event. This allows to attribute a Monte Carlo simulated event to a specific model (#1483). +Add GCTAEventList::set_mc_id_names() method to store Monte Carlo metadata into an +event list and to write this information into the FITS header. Furthermore, the +information is read back when loading an event list (#1483). + 21. COMPTEL interface --------------------- diff --git a/inst/cta/include/GCTAEventList.hpp b/inst/cta/include/GCTAEventList.hpp index d2b20d4..a0e4345 100644 --- a/inst/cta/include/GCTAEventList.hpp +++ b/inst/cta/include/GCTAEventList.hpp @@ -103,6 +103,8 @@ public: const bool& has_phase() const; const bool& has_detxy() const; const bool& has_mc_id() const; + void set_mc_id_names(const std::vector& ids, + const std::vector& names); protected: // Protected methods @@ -112,19 +114,23 @@ protected: virtual void set_energies(void) { return; } virtual void set_times(void) { return; } void read_events(const GFitsTable& table) const; + void read_mc_ids(const GFitsTable& table); void write_events(GFitsBinTable& table) const; void write_ds_keys(GFitsHDU& hdu, const std::string& gtiname = gammalib::extname_gti) const; + void write_mc_ids(GFitsHDU& hdu) const; int irf_cache_init(const std::string& name) const; int irf_cache_index(const std::string& name) const; // Event list meta data - GCTARoi m_roi; //!< Region of interest - int m_num_events; //!< Number of events - std::string m_gti_extname; //!< GTI extension name - bool m_has_phase; //!< Signal presence of phase - bool m_has_detxy; //!< Signal presence of detector coordinates - bool m_has_mc_id; //!< Signal presence of MC identifier + GCTARoi m_roi; //!< Region of interest + int m_num_events; //!< Number of events + std::string m_gti_extname; //!< GTI extension name + bool m_has_phase; //!< Signal presence of phase + bool m_has_detxy; //!< Signal presence of detector coordinates + bool m_has_mc_id; //!< Signal presence of MC identifier + std::vector m_mc_ids; //!< List of Monte Carlo IDs + std::vector m_mc_id_names; //!< List of model names // Event list data mutable std::vector m_events; //!< Events diff --git a/inst/cta/pyext/GCTAEventList.i b/inst/cta/pyext/GCTAEventList.i index 4f1f054..cbf4876 100644 --- a/inst/cta/pyext/GCTAEventList.i +++ b/inst/cta/pyext/GCTAEventList.i @@ -76,6 +76,8 @@ public: const bool& has_phase() const; const bool& has_detxy() const; const bool& has_mc_id() const; + void set_mc_id_names(const std::vector& ids, + const std::vector& names); }; diff --git a/inst/cta/src/GCTAEventAtom.cpp b/inst/cta/src/GCTAEventAtom.cpp index 5a4c776..987dfb1 100644 --- a/inst/cta/src/GCTAEventAtom.cpp +++ b/inst/cta/src/GCTAEventAtom.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAEventAtom.cpp - CTA event atom class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2015 by Jurgen Knodlseder * + * copyright (C) 2010-2017 by Jurgen Knodlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -209,6 +209,7 @@ void GCTAEventAtom::init_members(void) m_energy.clear(); m_index = 0; m_event_id = 0; + m_mc_id = 0; m_phase = 0.0; // Return @@ -229,6 +230,7 @@ void GCTAEventAtom::copy_members(const GCTAEventAtom& atom) m_energy = atom.m_energy; m_index = atom.m_index; m_event_id = atom.m_event_id; + m_mc_id = atom.m_mc_id; m_phase = atom.m_phase; // Return diff --git a/inst/cta/src/GCTAEventList.cpp b/inst/cta/src/GCTAEventList.cpp index 2429998..7ab3df6 100644 --- a/inst/cta/src/GCTAEventList.cpp +++ b/inst/cta/src/GCTAEventList.cpp @@ -49,6 +49,8 @@ #define G_OPERATOR "GCTAEventList::operator[](int&)" #define G_ROI "GCTAEventList::roi(GRoi&)" #define G_APPEND_COLUMN "GCTAEventList::append_column(GFitsTableCol&)" +#define G_SET_MC_ID_NAMES "GCTAEventList::set_mc_id_names(std::vector&,"\ + " std::vector&)" #define G_FETCH "GCTAEventList::fetch()" /* __ Macros _____________________________________________________________ */ @@ -442,6 +444,9 @@ void GCTAEventList::read(const GFits& fits) read_events(events); } + // Read Monte Carlo identifier keywords + read_mc_ids(events); + // Return return; } @@ -508,6 +513,9 @@ void GCTAEventList::write(GFits& fits, const std::string& evtname, // Write data selection keywords write_ds_keys(table, gtiname); + // Write Monte Carlo identifier keywords + write_mc_ids(table); + // Append event table to FITS file fits.append(table); @@ -792,9 +800,42 @@ void GCTAEventList::dispose(void) const /***********************************************************************//** + * @brief Set Monte Carlo identifiers and model names + * + * @param[in] ids Monte Carlo identifiers. + * @param[in] names Model name for each Monte Carlo identifier. + * + * Sets the Monte Carlo identifiers and their corresponding model names. + * The information will then be written into the header file. + ***************************************************************************/ +void GCTAEventList::set_mc_id_names(const std::vector& ids, + const std::vector& names) +{ + // Throw an exception if the vectors have an incompatible size + if (ids.size() != names.size()) { + std::string msg = "The number of Monte Carlo identifiers ("+ + gammalib::str(ids.size())+") does not correspond to " + "the number of model names ("+ + gammalib::str(names.size())+"). Please provide the " + "same number of model names and Monte Carlo " + "identifiers."; + throw GException::invalid_argument(G_SET_MC_ID_NAMES, msg); + } + + // Set vectors + m_mc_ids = ids; + m_mc_id_names = names; + + // Return + return; +} + + + +/***********************************************************************//** * @brief Print event list information * - * @param[in] chatter Chattiness (defaults to NORMAL). + * @param[in] chatter Chattiness. * @return String containing event list information. ***************************************************************************/ std::string GCTAEventList::print(const GChatter& chatter) const @@ -923,6 +964,8 @@ void GCTAEventList::init_members(void) m_has_phase = false; m_has_detxy = true; m_has_mc_id = false; + m_mc_ids.clear(); + m_mc_id_names.clear(); // Initialise cache m_irf_names.clear(); @@ -949,6 +992,8 @@ void GCTAEventList::copy_members(const GCTAEventList& list) m_has_phase = list.m_has_phase; m_has_detxy = list.m_has_detxy; m_has_mc_id = list.m_has_mc_id; + m_mc_ids = list.m_mc_ids; + m_mc_id_names = list.m_mc_id_names; // Copy cache m_irf_names = list.m_irf_names; @@ -1107,6 +1152,52 @@ void GCTAEventList::read_events(const GFitsTable& table) const /***********************************************************************//** + * @brief Read Monte Carlo identifier keywords from FITS HDU + * + * @param[in] table FITS table HDU. + * + * Reads the Monte Carlo identifier keywords for an event list from the FITS + * table HDU. + ***************************************************************************/ +void GCTAEventList::read_mc_ids(const GFitsTable& table) +{ + // Clear Monte Carlo identifiers and names + m_mc_ids.clear(); + m_mc_id_names.clear(); + + // Continue only if the header contains a "MCIDS" keyword + if (table.has_card("NMCIDS")) { + + // Get number of Monte Carlo identifiers + int nids = table.integer("NMCIDS"); + + // Loop over Monte Carlo identifiers + for (int i = 0; i < nids; ++i) { + + // Set keyword names + char keyword_id[10]; + char keyword_name[10]; + sprintf(keyword_id, "MID%5.5d", i+1); + sprintf(keyword_name, "MMN%5.5d", i+1); + + // Get header keywords + int id = table.integer("keyword_id"); + std::string name = table.string("keyword_name"); + + // Put identifier and name in list + m_mc_ids.push_back(id); + m_mc_id_names.push_back(name); + + } // endfor: looped over Monte Carlo identifiers + + } // endif: there were Monte Carlo identifiers + + // Return + return; +} + + +/***********************************************************************//** * @brief Write CTA events into FITS table * * @param[in] hdu FITS binary table. @@ -1248,16 +1339,17 @@ void GCTAEventList::write_events(GFitsBinTable& hdu) const * Writes the data sub-space keywords for an event list into the FITS HDU. * The following keywords will be written: * - * DSTYP1 = "TIME" / Data sub-space type - * DSUNI1 = "s" / Data sub-space unit - * DSVAL1 = "TABLE" / Data sub-space value - * DSREF1 = ":[extname]" / Data sub-space reference - * DSTYP2 = "ENERGY" / Data sub-space type - * DSUNI2 = "TeV" / Data sub-space unit - * DSVAL2 = "[emin]:[emax]" / Data sub-space value - * DSTYP3 = "POS(RA,DEC)" / Data sub-space type - * DSUNI3 = "deg" / Data sub-space unit - * DSVAL3 = "CIRCLE([ra],[dec],[rad])" / Data sub-space value + * DSTYP1 = "TIME" / Data sub-space type + * DSUNI1 = "s" / Data sub-space unit + * DSVAL1 = "TABLE" / Data sub-space value + * DSREF1 = ":[extname]" / Data sub-space reference + * DSTYP2 = "ENERGY" / Data sub-space type + * DSUNI2 = "TeV" / Data sub-space unit + * DSVAL2 = "[emin]:[emax]" / Data sub-space value + * DSTYP3 = "POS(RA,DEC)" / Data sub-space type + * DSUNI3 = "deg" / Data sub-space unit + * DSVAL3 = "CIRCLE([ra],[dec],[rad])" / Data sub-space value + * NDSKEYS = 3 / Number of data sub-space keys * * where * @@ -1326,6 +1418,59 @@ void GCTAEventList::write_ds_keys(GFitsHDU& hdu, const std::string& gtiname) con /***********************************************************************//** + * @brief Write Monte Carlo identifier keywords into FITS HDU + * + * @param[in] hdu FITS HDU. + * + * Writes the Monte Carlo identifier keywords for an event list into the FITS + * HDU. The following keywords will be written: + * + * NMCIDS Number of Monte Carlo identifiers + * MID00001 First Monte Carlo identifier + * MMN00001 Model name for first Monte Carlo identifier + * MID00002 Second Monte Carlo identifier + * MMN00002 Model name for second Monte Carlo identifier + * ... + ***************************************************************************/ +void GCTAEventList::write_mc_ids(GFitsHDU& hdu) const +{ + // Set number of Monte Carlo identifiers + int nids = m_mc_ids.size(); + + // Continue only if there are Monte Carlo identifiers + if (nids > 0) { + + // Set number of Monte Carlo identifiers + hdu.card("NMCIDS", nids, "Number of Monte Carlo identifiers"); + + // Loop over Monte Carlo identifiers + for (int i = 0; i < nids; ++i) { + + // Set keyword names + char keyword_id[10]; + char keyword_name[10]; + sprintf(keyword_id, "MID%5.5d", i+1); + sprintf(keyword_name, "MMN%5.5d", i+1); + + // Set comments + std::string comment_id = "Monte Carlo identifier for model " + + gammalib::str(i+1); + std::string comment_name = "Name of model " + gammalib::str(i+1); + + // Write header keywords + hdu.card(std::string(keyword_id), m_mc_ids[i], comment_id); + hdu.card(std::string(keyword_name), m_mc_id_names[i], comment_name); + + } // endfor: looped over Monte Carlo identifiers + + } // endif: there were Monte Carlo identifiers + + // Return + return; +} + + +/***********************************************************************//** * @brief Initialize IRF cache for a given model * * @param[in] name Model name. -- libgit2 0.21.4