Peter Deiml / gammalib

Commit 8f3f4932f78b27f75a684c5d8bc3eb8aa37ca6aa

Authored by Jürgen Knödlseder
1 parent b0375e86

Add support for Monte Carlo identifier metadata (#1483)

Inline Side-by-side
Showing 6 changed files with 181 additions and 20 deletions
ChangeLog
  View file @8f3f493
1   -2017-04-13
  1 +2017-04-14
2 2  
3 3 * Version 1.3.0 released
4 4 =========================
5 5  
  6 + Add support for Monte Carlo identifier metadata (#1483)
6 7 Add GEventList::has_mc_id() and GEventAtom::mc_id() methods (#1483)
7 8 Add GEventList::has_phase() and GEventList::has_detxy() methods (#1982)
8 9 Add GEventList::append_column() method (#1991)
... ...
1 1 New Features and Important Changes in GammaLib 1.3.0.dev1
2 2  
3   -13 April 2017
  3 +14 April 2017
4 4  
5 5  
6 6 1. Introduction
... ... @@ -54,6 +54,7 @@ The following methods have been added:
54 54 - GCTAEventList::has_detxy()
55 55 - GCTAEventList::has_mc_id(bool&)
56 56 - GCTAEventList::has_mc_id()
  57 +- GCTAEventList::set_mc_id_names(std::vector<int>&, std::vector<std::string>&)
57 58 - GCTAEventAtom::mc_id(int&)
58 59 - GCTAEventAtom::mc_id()
59 60  
... ... @@ -233,6 +234,10 @@ Add GCTAEventList::has_mc_id() and GCTAEventAtom::mc_id() getter and setter meth
233 234 to set the Monte Carlo identifier of a given event. This allows to attribute a
234 235 Monte Carlo simulated event to a specific model (#1483).
235 236  
  237 +Add GCTAEventList::set_mc_id_names() method to store Monte Carlo metadata into an
  238 +event list and to write this information into the FITS header. Furthermore, the
  239 +information is read back when loading an event list (#1483).
  240 +
236 241  
237 242 21. COMPTEL interface
238 243 ---------------------
... ...
inst/cta/include/GCTAEventList.hpp
  View file @8f3f493
... ... @@ -103,6 +103,8 @@ public:
103 103 const bool& has_phase() const;
104 104 const bool& has_detxy() const;
105 105 const bool& has_mc_id() const;
  106 + void set_mc_id_names(const std::vector<int>& ids,
  107 + const std::vector<std::string>& names);
106 108  
107 109 protected:
108 110 // Protected methods
... ... @@ -112,19 +114,23 @@ protected:
112 114 virtual void set_energies(void) { return; }
113 115 virtual void set_times(void) { return; }
114 116 void read_events(const GFitsTable& table) const;
  117 + void read_mc_ids(const GFitsTable& table);
115 118 void write_events(GFitsBinTable& table) const;
116 119 void write_ds_keys(GFitsHDU& hdu,
117 120 const std::string& gtiname = gammalib::extname_gti) const;
  121 + void write_mc_ids(GFitsHDU& hdu) const;
118 122 int irf_cache_init(const std::string& name) const;
119 123 int irf_cache_index(const std::string& name) const;
120 124  
121 125 // Event list meta data
122   - GCTARoi m_roi; //!< Region of interest
123   - int m_num_events; //!< Number of events
124   - std::string m_gti_extname; //!< GTI extension name
125   - bool m_has_phase; //!< Signal presence of phase
126   - bool m_has_detxy; //!< Signal presence of detector coordinates
127   - bool m_has_mc_id; //!< Signal presence of MC identifier
  126 + GCTARoi m_roi; //!< Region of interest
  127 + int m_num_events; //!< Number of events
  128 + std::string m_gti_extname; //!< GTI extension name
  129 + bool m_has_phase; //!< Signal presence of phase
  130 + bool m_has_detxy; //!< Signal presence of detector coordinates
  131 + bool m_has_mc_id; //!< Signal presence of MC identifier
  132 + std::vector<int> m_mc_ids; //!< List of Monte Carlo IDs
  133 + std::vector<std::string> m_mc_id_names; //!< List of model names
128 134  
129 135 // Event list data
130 136 mutable std::vector<GCTAEventAtom> m_events; //!< Events
... ...
inst/cta/pyext/GCTAEventList.i
  View file @8f3f493
... ... @@ -76,6 +76,8 @@ public:
76 76 const bool& has_phase() const;
77 77 const bool& has_detxy() const;
78 78 const bool& has_mc_id() const;
  79 + void set_mc_id_names(const std::vector<int>& ids,
  80 + const std::vector<std::string>& names);
79 81 };
80 82  
81 83  
... ...
inst/cta/src/GCTAEventAtom.cpp
  View file @8f3f493
1 1 /***************************************************************************
2 2 * GCTAEventAtom.cpp - CTA event atom class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2015 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2017 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -209,6 +209,7 @@ void GCTAEventAtom::init_members(void)
209 209 m_energy.clear();
210 210 m_index = 0;
211 211 m_event_id = 0;
  212 + m_mc_id = 0;
212 213 m_phase = 0.0;
213 214  
214 215 // Return
... ... @@ -229,6 +230,7 @@ void GCTAEventAtom::copy_members(const GCTAEventAtom&amp; atom)
229 230 m_energy = atom.m_energy;
230 231 m_index = atom.m_index;
231 232 m_event_id = atom.m_event_id;
  233 + m_mc_id = atom.m_mc_id;
232 234 m_phase = atom.m_phase;
233 235  
234 236 // Return
... ...
inst/cta/src/GCTAEventList.cpp
  View file @8f3f493
... ... @@ -49,6 +49,8 @@
49 49 #define G_OPERATOR "GCTAEventList::operator[](int&)"
50 50 #define G_ROI "GCTAEventList::roi(GRoi&)"
51 51 #define G_APPEND_COLUMN "GCTAEventList::append_column(GFitsTableCol&)"
  52 +#define G_SET_MC_ID_NAMES "GCTAEventList::set_mc_id_names(std::vector<int>&,"\
  53 + " std::vector<std::string>&)"
52 54 #define G_FETCH "GCTAEventList::fetch()"
53 55  
54 56 /* __ Macros _____________________________________________________________ */
... ... @@ -442,6 +444,9 @@ void GCTAEventList::read(const GFits&amp; fits)
442 444 read_events(events);
443 445 }
444 446  
  447 + // Read Monte Carlo identifier keywords
  448 + read_mc_ids(events);
  449 +
445 450 // Return
446 451 return;
447 452 }
... ... @@ -508,6 +513,9 @@ void GCTAEventList::write(GFits&amp; fits, const std::string&amp; evtname,
508 513 // Write data selection keywords
509 514 write_ds_keys(table, gtiname);
510 515  
  516 + // Write Monte Carlo identifier keywords
  517 + write_mc_ids(table);
  518 +
511 519 // Append event table to FITS file
512 520 fits.append(table);
513 521  
... ... @@ -792,9 +800,42 @@ void GCTAEventList::dispose(void) const
792 800  
793 801  
794 802 /***********************************************************************//**
  803 + * @brief Set Monte Carlo identifiers and model names
  804 + *
  805 + * @param[in] ids Monte Carlo identifiers.
  806 + * @param[in] names Model name for each Monte Carlo identifier.
  807 + *
  808 + * Sets the Monte Carlo identifiers and their corresponding model names.
  809 + * The information will then be written into the header file.
  810 + ***************************************************************************/
  811 +void GCTAEventList::set_mc_id_names(const std::vector<int>& ids,
  812 + const std::vector<std::string>& names)
  813 +{
  814 + // Throw an exception if the vectors have an incompatible size
  815 + if (ids.size() != names.size()) {
  816 + std::string msg = "The number of Monte Carlo identifiers ("+
  817 + gammalib::str(ids.size())+") does not correspond to "
  818 + "the number of model names ("+
  819 + gammalib::str(names.size())+"). Please provide the "
  820 + "same number of model names and Monte Carlo "
  821 + "identifiers.";
  822 + throw GException::invalid_argument(G_SET_MC_ID_NAMES, msg);
  823 + }
  824 +
  825 + // Set vectors
  826 + m_mc_ids = ids;
  827 + m_mc_id_names = names;
  828 +
  829 + // Return
  830 + return;
  831 +}
  832 +
  833 +
  834 +
  835 +/***********************************************************************//**
795 836 * @brief Print event list information
796 837 *
797   - * @param[in] chatter Chattiness (defaults to NORMAL).
  838 + * @param[in] chatter Chattiness.
798 839 * @return String containing event list information.
799 840 ***************************************************************************/
800 841 std::string GCTAEventList::print(const GChatter& chatter) const
... ... @@ -923,6 +964,8 @@ void GCTAEventList::init_members(void)
923 964 m_has_phase = false;
924 965 m_has_detxy = true;
925 966 m_has_mc_id = false;
  967 + m_mc_ids.clear();
  968 + m_mc_id_names.clear();
926 969  
927 970 // Initialise cache
928 971 m_irf_names.clear();
... ... @@ -949,6 +992,8 @@ void GCTAEventList::copy_members(const GCTAEventList&amp; list)
949 992 m_has_phase = list.m_has_phase;
950 993 m_has_detxy = list.m_has_detxy;
951 994 m_has_mc_id = list.m_has_mc_id;
  995 + m_mc_ids = list.m_mc_ids;
  996 + m_mc_id_names = list.m_mc_id_names;
952 997  
953 998 // Copy cache
954 999 m_irf_names = list.m_irf_names;
... ... @@ -1107,6 +1152,52 @@ void GCTAEventList::read_events(const GFitsTable&amp; table) const
1107 1152  
1108 1153  
1109 1154 /***********************************************************************//**
  1155 + * @brief Read Monte Carlo identifier keywords from FITS HDU
  1156 + *
  1157 + * @param[in] table FITS table HDU.
  1158 + *
  1159 + * Reads the Monte Carlo identifier keywords for an event list from the FITS
  1160 + * table HDU.
  1161 + ***************************************************************************/
  1162 +void GCTAEventList::read_mc_ids(const GFitsTable& table)
  1163 +{
  1164 + // Clear Monte Carlo identifiers and names
  1165 + m_mc_ids.clear();
  1166 + m_mc_id_names.clear();
  1167 +
  1168 + // Continue only if the header contains a "MCIDS" keyword
  1169 + if (table.has_card("NMCIDS")) {
  1170 +
  1171 + // Get number of Monte Carlo identifiers
  1172 + int nids = table.integer("NMCIDS");
  1173 +
  1174 + // Loop over Monte Carlo identifiers
  1175 + for (int i = 0; i < nids; ++i) {
  1176 +
  1177 + // Set keyword names
  1178 + char keyword_id[10];
  1179 + char keyword_name[10];
  1180 + sprintf(keyword_id, "MID%5.5d", i+1);
  1181 + sprintf(keyword_name, "MMN%5.5d", i+1);
  1182 +
  1183 + // Get header keywords
  1184 + int id = table.integer("keyword_id");
  1185 + std::string name = table.string("keyword_name");
  1186 +
  1187 + // Put identifier and name in list
  1188 + m_mc_ids.push_back(id);
  1189 + m_mc_id_names.push_back(name);
  1190 +
  1191 + } // endfor: looped over Monte Carlo identifiers
  1192 +
  1193 + } // endif: there were Monte Carlo identifiers
  1194 +
  1195 + // Return
  1196 + return;
  1197 +}
  1198 +
  1199 +
  1200 +/***********************************************************************//**
1110 1201 * @brief Write CTA events into FITS table
1111 1202 *
1112 1203 * @param[in] hdu FITS binary table.
... ... @@ -1248,16 +1339,17 @@ void GCTAEventList::write_events(GFitsBinTable&amp; hdu) const
1248 1339 * Writes the data sub-space keywords for an event list into the FITS HDU.
1249 1340 * The following keywords will be written:
1250 1341 *
1251   - * DSTYP1 = "TIME" / Data sub-space type
1252   - * DSUNI1 = "s" / Data sub-space unit
1253   - * DSVAL1 = "TABLE" / Data sub-space value
1254   - * DSREF1 = ":[extname]" / Data sub-space reference
1255   - * DSTYP2 = "ENERGY" / Data sub-space type
1256   - * DSUNI2 = "TeV" / Data sub-space unit
1257   - * DSVAL2 = "[emin]:[emax]" / Data sub-space value
1258   - * DSTYP3 = "POS(RA,DEC)" / Data sub-space type
1259   - * DSUNI3 = "deg" / Data sub-space unit
1260   - * DSVAL3 = "CIRCLE([ra],[dec],[rad])" / Data sub-space value
  1342 + * DSTYP1 = "TIME" / Data sub-space type
  1343 + * DSUNI1 = "s" / Data sub-space unit
  1344 + * DSVAL1 = "TABLE" / Data sub-space value
  1345 + * DSREF1 = ":[extname]" / Data sub-space reference
  1346 + * DSTYP2 = "ENERGY" / Data sub-space type
  1347 + * DSUNI2 = "TeV" / Data sub-space unit
  1348 + * DSVAL2 = "[emin]:[emax]" / Data sub-space value
  1349 + * DSTYP3 = "POS(RA,DEC)" / Data sub-space type
  1350 + * DSUNI3 = "deg" / Data sub-space unit
  1351 + * DSVAL3 = "CIRCLE([ra],[dec],[rad])" / Data sub-space value
  1352 + * NDSKEYS = 3 / Number of data sub-space keys
1261 1353 *
1262 1354 * where
1263 1355 *
... ... @@ -1326,6 +1418,59 @@ void GCTAEventList::write_ds_keys(GFitsHDU&amp; hdu, const std::string&amp; gtiname) con
1326 1418  
1327 1419  
1328 1420 /***********************************************************************//**
  1421 + * @brief Write Monte Carlo identifier keywords into FITS HDU
  1422 + *
  1423 + * @param[in] hdu FITS HDU.
  1424 + *
  1425 + * Writes the Monte Carlo identifier keywords for an event list into the FITS
  1426 + * HDU. The following keywords will be written:
  1427 + *
  1428 + * NMCIDS Number of Monte Carlo identifiers
  1429 + * MID00001 First Monte Carlo identifier
  1430 + * MMN00001 Model name for first Monte Carlo identifier
  1431 + * MID00002 Second Monte Carlo identifier
  1432 + * MMN00002 Model name for second Monte Carlo identifier
  1433 + * ...
  1434 + ***************************************************************************/
  1435 +void GCTAEventList::write_mc_ids(GFitsHDU& hdu) const
  1436 +{
  1437 + // Set number of Monte Carlo identifiers
  1438 + int nids = m_mc_ids.size();
  1439 +
  1440 + // Continue only if there are Monte Carlo identifiers
  1441 + if (nids > 0) {
  1442 +
  1443 + // Set number of Monte Carlo identifiers
  1444 + hdu.card("NMCIDS", nids, "Number of Monte Carlo identifiers");
  1445 +
  1446 + // Loop over Monte Carlo identifiers
  1447 + for (int i = 0; i < nids; ++i) {
  1448 +
  1449 + // Set keyword names
  1450 + char keyword_id[10];
  1451 + char keyword_name[10];
  1452 + sprintf(keyword_id, "MID%5.5d", i+1);
  1453 + sprintf(keyword_name, "MMN%5.5d", i+1);
  1454 +
  1455 + // Set comments
  1456 + std::string comment_id = "Monte Carlo identifier for model " +
  1457 + gammalib::str(i+1);
  1458 + std::string comment_name = "Name of model " + gammalib::str(i+1);
  1459 +
  1460 + // Write header keywords
  1461 + hdu.card(std::string(keyword_id), m_mc_ids[i], comment_id);
  1462 + hdu.card(std::string(keyword_name), m_mc_id_names[i], comment_name);
  1463 +
  1464 + } // endfor: looped over Monte Carlo identifiers
  1465 +
  1466 + } // endif: there were Monte Carlo identifiers
  1467 +
  1468 + // Return
  1469 + return;
  1470 +}
  1471 +
  1472 +
  1473 +/***********************************************************************//**
1329 1474 * @brief Initialize IRF cache for a given model
1330 1475 *
1331 1476 * @param[in] name Model name.
... ...