Francois Brun / gammalib

Commit d1842e6f4bea3ddac7735cd948b3efa2e3ff7aa5

Authored by Jürgen Knödlseder
2 parents 7f2ae9e0 cbde4004

Merge branch '1600-GCTAEventList-flexible-reading-writing' into integration

Inline Side-by-side
Showing 7 changed files with 158 additions and 425 deletions
ChangeLog
  View file @d1842e6
1 -2015-12-02 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu> 1 +2015-12-16 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2
3 * Version 1.1.0 released 3 * Version 1.1.0 released
4 ========================= 4 =========================
5 5
  6 + Read only mandatary columns in GCTAEventList (#1600)
6 Add GModel::has_ts() method (#1578) 7 Add GModel::has_ts() method (#1578)
7 8
8 9
@@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0 @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0
2 2
3 Author(s) : Juergen Knoedlseder 3 Author(s) : Juergen Knoedlseder
4 4
5 -2 December 2015 5 +16 December 2015
6 6
7 7
8 1. Introduction 8 1. Introduction
@@ -109,9 +109,12 @@ None @@ -109,9 +109,12 @@ None
109 109
110 20. CTA interface 110 20. CTA interface
111 ----------------- 111 -----------------
112 -The IRF loaders now check the existence of the file before loading, and now throw 112 +The IRF loaders now check the existence of the file before loading, and now throws
113 appropriate error messages in case that the files do not exist (#1586). 113 appropriate error messages in case that the files do not exist (#1586).
114 114
  115 +Only manadatory columns are read by GCTAEventsList, while any other columns will only
  116 +be stored for back writing in a FITS file.
  117 +
115 118
116 21. COMPTEL interface 119 21. COMPTEL interface
117 --------------------- 120 ---------------------
inst/cta/include/GCTAEventAtom.hpp
  View file @d1842e6
@@ -72,10 +72,8 @@ public: @@ -72,10 +72,8 @@ public:
72 // Other methods 72 // Other methods
73 const int& index(void) const; 73 const int& index(void) const;
74 const unsigned long& event_id(void) const; 74 const unsigned long& event_id(void) const;
75 - const unsigned long& obs_id(void) const;  
76 const float& phase(void) const; 75 const float& phase(void) const;
77 void event_id(const unsigned long& id); 76 void event_id(const unsigned long& id);
78 - void obs_id(const unsigned long& id);  
79 void phase(const float& phase); 77 void phase(const float& phase);
80 78
81 protected: 79 protected:
@@ -90,24 +88,6 @@ protected: @@ -90,24 +88,6 @@ protected:
90 GEnergy m_energy; //!< Event energy 88 GEnergy m_energy; //!< Event energy
91 GTime m_time; //!< Event time 89 GTime m_time; //!< Event time
92 unsigned long m_event_id; //!< Event identifier 90 unsigned long m_event_id; //!< Event identifier
93 - unsigned long m_obs_id; //!< Observation identifier  
94 - int m_multip; //!< Multiplicity  
95 - char m_telmask; //!< Telescope mask  
96 - float m_dir_err; //!< Error on event direction  
97 - float m_alt; //!< Event altitude  
98 - float m_az; //!< Event azimuth  
99 - float m_corex; //!< Position on ground (m)  
100 - float m_corey; //!< Position on ground (m)  
101 - float m_core_err; //!< Error on core reconstruction  
102 - float m_xmax; //!< Position of shower max (g/cm2)  
103 - float m_xmax_err; //!< Error on shower max (g/cm2)  
104 - float m_shwidth; //!< Shower width (m)  
105 - float m_shlength; //!< Shower length (m)  
106 - float m_energy_err; //!< Error on event energy (MeV)  
107 - float m_hil_msw; //!< Hillas MSW  
108 - float m_hil_msw_err; //!< Hillas MSW error  
109 - float m_hil_msl; //!< Hillas MSL  
110 - float m_hil_msl_err; //!< Hillas MSL error  
111 float m_phase; //!< Optional phase 91 float m_phase; //!< Optional phase
112 }; 92 };
113 93
@@ -226,18 +206,6 @@ const unsigned long&amp; GCTAEventAtom::event_id(void) const @@ -226,18 +206,6 @@ const unsigned long&amp; GCTAEventAtom::event_id(void) const
226 206
227 207
228 /***********************************************************************//** 208 /***********************************************************************//**
229 - * @brief Return observation identifier  
230 - *  
231 - * @return Observation identifier.  
232 - ***************************************************************************/  
233 -inline  
234 -const unsigned long& GCTAEventAtom::obs_id(void) const  
235 -{  
236 - return (m_obs_id);  
237 -}  
238 -  
239 -  
240 -/***********************************************************************//**  
241 * @brief Return event phase 209 * @brief Return event phase
242 * 210 *
243 * @return Event phase. 211 * @return Event phase.
@@ -263,19 +231,6 @@ void GCTAEventAtom::event_id(const unsigned long&amp; id) @@ -263,19 +231,6 @@ void GCTAEventAtom::event_id(const unsigned long&amp; id)
263 231
264 232
265 /***********************************************************************//** 233 /***********************************************************************//**
266 - * @brief Set observation identifier  
267 - *  
268 - * @param[in] id Observation identifier.  
269 - ***************************************************************************/  
270 -inline  
271 -void GCTAEventAtom::obs_id(const unsigned long& id)  
272 -{  
273 - m_obs_id = id;  
274 - return;  
275 -}  
276 -  
277 -  
278 -/***********************************************************************//**  
279 * @brief Set event phase 234 * @brief Set event phase
280 * 235 *
281 * @param[in] phase Event phase. 236 * @param[in] phase Event phase.
inst/cta/include/GCTAEventList.hpp
  View file @d1842e6
@@ -88,19 +88,18 @@ protected: @@ -88,19 +88,18 @@ protected:
88 void free_members(void); 88 void free_members(void);
89 virtual void set_energies(void) { return; } 89 virtual void set_energies(void) { return; }
90 virtual void set_times(void) { return; } 90 virtual void set_times(void) { return; }
91 - void read_events(const GFitsTable& hdu);  
92 - void read_events_v0(const GFitsTable& hdu);  
93 - void read_events_v1(const GFitsTable& hdu);  
94 - void read_events_hillas(const GFitsTable& hdu);  
95 - void write_events(GFitsBinTable& hdu) const; 91 + void read_events(const GFitsTable& table);
  92 + void write_events(GFitsBinTable& table) const;
96 void write_ds_keys(GFitsHDU& hdu) const; 93 void write_ds_keys(GFitsHDU& hdu) const;
97 int irf_cache_init(const std::string& name) const; 94 int irf_cache_init(const std::string& name) const;
98 int irf_cache_index(const std::string& name) const; 95 int irf_cache_index(const std::string& name) const;
99 96
100 // Protected members 97 // Protected members
101 - GCTARoi m_roi; //!< Region of interest  
102 - std::vector<GCTAEventAtom> m_events; //!< Events  
103 - bool m_has_phase; //!< Signal presence of phase 98 + GCTARoi m_roi; //!< Region of interest
  99 + std::vector<GCTAEventAtom> m_events; //!< Events
  100 + std::vector<GFitsTableCol*> m_columns; //!< Pointers to optional columns
  101 + bool m_has_phase; //!< Signal presence of phase
  102 + bool m_has_detxy; //!< Signal presence of detector coordinates
104 103
105 // IRF cache for diffuse models 104 // IRF cache for diffuse models
106 mutable std::vector<std::string> m_irf_names; //!< Model names 105 mutable std::vector<std::string> m_irf_names; //!< Model names
inst/cta/pyext/GCTAEventAtom.i
  View file @d1842e6
@@ -55,10 +55,8 @@ public: @@ -55,10 +55,8 @@ public:
55 // Other methods 55 // Other methods
56 const int& index(void) const; 56 const int& index(void) const;
57 const unsigned long& event_id(void) const; 57 const unsigned long& event_id(void) const;
58 - const unsigned long& obs_id(void) const;  
59 const float& phase(void) const; 58 const float& phase(void) const;
60 void event_id(const unsigned long& id); 59 void event_id(const unsigned long& id);
61 - void obs_id(const unsigned long& id);  
62 void phase(const float& phase); 60 void phase(const float& phase);
63 }; 61 };
64 62
inst/cta/src/GCTAEventAtom.cpp
  View file @d1842e6
@@ -207,27 +207,9 @@ void GCTAEventAtom::init_members(void) @@ -207,27 +207,9 @@ void GCTAEventAtom::init_members(void)
207 m_dir.clear(); 207 m_dir.clear();
208 m_time.clear(); 208 m_time.clear();
209 m_energy.clear(); 209 m_energy.clear();
210 - m_index = 0;  
211 - m_event_id = 0;  
212 - m_obs_id = 0;  
213 - m_multip = 0;  
214 - m_telmask = 0;  
215 - m_dir_err = 0.0;  
216 - m_alt = 0.0;  
217 - m_az = 0.0;  
218 - m_corex = 0.0;  
219 - m_corey = 0.0;  
220 - m_core_err = 0.0;  
221 - m_xmax = 0.0;  
222 - m_xmax_err = 0.0;  
223 - m_shwidth = 0.0;  
224 - m_shlength = 0.0;  
225 - m_energy_err = 0.0;  
226 - m_hil_msw = 0.0;  
227 - m_hil_msw_err = 0.0;  
228 - m_hil_msl = 0.0;  
229 - m_hil_msl_err = 0.0;  
230 - m_phase = 0.0; 210 + m_index = 0;
  211 + m_event_id = 0;
  212 + m_phase = 0.0;
231 213
232 // Return 214 // Return
233 return; 215 return;
@@ -242,30 +224,12 @@ void GCTAEventAtom::init_members(void) @@ -242,30 +224,12 @@ void GCTAEventAtom::init_members(void)
242 void GCTAEventAtom::copy_members(const GCTAEventAtom& atom) 224 void GCTAEventAtom::copy_members(const GCTAEventAtom& atom)
243 { 225 {
244 // Copy members 226 // Copy members
245 - m_dir = atom.m_dir;  
246 - m_time = atom.m_time;  
247 - m_energy = atom.m_energy;  
248 - m_index = atom.m_index;  
249 - m_event_id = atom.m_event_id;  
250 - m_obs_id = atom.m_obs_id;  
251 - m_multip = atom.m_multip;  
252 - m_telmask = atom.m_telmask;  
253 - m_dir_err = atom.m_dir_err;  
254 - m_alt = atom.m_alt;  
255 - m_az = atom.m_az;  
256 - m_corex = atom.m_corex;  
257 - m_corey = atom.m_corey;  
258 - m_core_err = atom.m_core_err;  
259 - m_xmax = atom.m_xmax;  
260 - m_xmax_err = atom.m_xmax_err;  
261 - m_shwidth = atom.m_shwidth;  
262 - m_shlength = atom.m_shlength;  
263 - m_energy_err = atom.m_energy_err;  
264 - m_hil_msw = atom.m_hil_msw;  
265 - m_hil_msw_err = atom.m_hil_msw_err;  
266 - m_hil_msl = atom.m_hil_msl;  
267 - m_hil_msl_err = atom.m_hil_msl_err;  
268 - m_phase = atom.m_phase; 227 + m_dir = atom.m_dir;
  228 + m_time = atom.m_time;
  229 + m_energy = atom.m_energy;
  230 + m_index = atom.m_index;
  231 + m_event_id = atom.m_event_id;
  232 + m_phase = atom.m_phase;
269 233
270 // Return 234 // Return
271 return; 235 return;
inst/cta/src/GCTAEventList.cpp
  View file @d1842e6
@@ -572,7 +572,9 @@ void GCTAEventList::init_members(void) @@ -572,7 +572,9 @@ void GCTAEventList::init_members(void)
572 // Initialise members 572 // Initialise members
573 m_roi.clear(); 573 m_roi.clear();
574 m_events.clear(); 574 m_events.clear();
  575 + m_columns.clear();
575 m_has_phase = false; 576 m_has_phase = false;
  577 + m_has_detxy = true;
576 578
577 // Initialise cache 579 // Initialise cache
578 m_irf_names.clear(); 580 m_irf_names.clear();
@@ -594,11 +596,18 @@ void GCTAEventList::copy_members(const GCTAEventList&amp; list) @@ -594,11 +596,18 @@ void GCTAEventList::copy_members(const GCTAEventList&amp; list)
594 m_roi = list.m_roi; 596 m_roi = list.m_roi;
595 m_events = list.m_events; 597 m_events = list.m_events;
596 m_has_phase = list.m_has_phase; 598 m_has_phase = list.m_has_phase;
  599 + m_has_detxy = list.m_has_detxy;
597 600
598 // Copy cache 601 // Copy cache
599 m_irf_names = list.m_irf_names; 602 m_irf_names = list.m_irf_names;
600 m_irf_values = list.m_irf_values; 603 m_irf_values = list.m_irf_values;
601 604
  605 + // Copy column pointers
  606 + m_columns.clear();
  607 + for (int i = 0; i < list.m_columns.size(); ++i) {
  608 + m_columns.push_back((list.m_columns[i]->clone()));
  609 + }
  610 +
602 // Return 611 // Return
603 return; 612 return;
604 } 613 }
@@ -609,6 +618,12 @@ void GCTAEventList::copy_members(const GCTAEventList&amp; list) @@ -609,6 +618,12 @@ void GCTAEventList::copy_members(const GCTAEventList&amp; list)
609 ***************************************************************************/ 618 ***************************************************************************/
610 void GCTAEventList::free_members(void) 619 void GCTAEventList::free_members(void)
611 { 620 {
  621 + // Free columns
  622 + for (int i = 0; i < m_columns.size(); ++i) {
  623 + delete m_columns[i];
  624 + m_columns[i] = NULL;
  625 + }
  626 +
612 // Return 627 // Return
613 return; 628 return;
614 } 629 }
@@ -619,52 +634,26 @@ void GCTAEventList::free_members(void) @@ -619,52 +634,26 @@ void GCTAEventList::free_members(void)
619 * 634 *
620 * @param[in] table FITS table. 635 * @param[in] table FITS table.
621 * 636 *
622 - * This method reads the CTA event list from a FITS table HDU into memory.  
623 - * Depending on the columns existing in the file, it either selects v0 or  
624 - * v1 of the event list reader.  
625 - ***************************************************************************/  
626 -void GCTAEventList::read_events(const GFitsTable& table)  
627 -{  
628 - // Clear existing events  
629 - m_events.clear();  
630 -  
631 - // Extract number of events in FITS file  
632 - int num = table.integer("NAXIS2");  
633 -  
634 - // Continue only if there are events  
635 - if (num > 0) {  
636 -  
637 - // Read events for v1  
638 - if (table.contains("SHWIDTH") && table.contains("SHLENGTH")) {  
639 - read_events_v1(table);  
640 - }  
641 -  
642 - // ... otherwise read events for v0  
643 - else {  
644 - read_events_v0(table);  
645 - }  
646 -  
647 - // Read (optional) Hillas parameters  
648 - read_events_hillas(table);  
649 -  
650 - } // endif: there were events  
651 -  
652 - // Return  
653 - return;  
654 -}  
655 -  
656 -  
657 -/***********************************************************************//**  
658 - * @brief Read CTA events from FITS table (version 0) 637 + * This method reads the CTA event list from a FITS table into memory.
  638 + * The following columns are mandatory in the FITS table:
659 * 639 *
660 - * @param[in] table FITS table. 640 + * - EVENT_ID
  641 + * - TIME
  642 + * - RA
  643 + * - DEC
  644 + * - ENERGY
  645 + *
  646 + * Optionally, the
661 * 647 *
662 - * This method reads the CTA event list from a FITS table HDU into memory.  
663 - * It is a minimal event reader that is compliant with the initial data  
664 - * format distributed by Karl Kosack. Information that is not present in  
665 - * that format is set to 0. 648 + * - DETX
  649 + * - DETY
  650 + * - PHASE
  651 + *
  652 + * columns are read. Any other columns present in the FITS table will also be
  653 + * read, and written into the FITS table when the write_events() method is
  654 + * called.
666 ***************************************************************************/ 655 ***************************************************************************/
667 -void GCTAEventList::read_events_v0(const GFitsTable& table) 656 +void GCTAEventList::read_events(const GFitsTable& table)
668 { 657 {
669 // Clear existing events 658 // Clear existing events
670 m_events.clear(); 659 m_events.clear();
@@ -679,120 +668,25 @@ void GCTAEventList::read_events_v0(const GFitsTable&amp; table) @@ -679,120 +668,25 @@ void GCTAEventList::read_events_v0(const GFitsTable&amp; table)
679 m_events.reserve(num); 668 m_events.reserve(num);
680 669
681 // Get column pointers 670 // Get column pointers
682 - const GFitsTableCol* ptr_eid = table["EVENT_ID"];  
683 - const GFitsTableCol* ptr_time = table["TIME"];  
684 - const GFitsTableCol* ptr_multip = table["MULTIP"];  
685 - const GFitsTableCol* ptr_ra = table["RA"];  
686 - const GFitsTableCol* ptr_dec = table["DEC"];  
687 - const GFitsTableCol* ptr_dir_err = table["DIR_ERR"];  
688 - const GFitsTableCol* ptr_detx = table["DETX"];  
689 - const GFitsTableCol* ptr_dety = table["DETY"];  
690 - const GFitsTableCol* ptr_alt = table["ALT"];  
691 - const GFitsTableCol* ptr_az = table["AZ"];  
692 - const GFitsTableCol* ptr_corex = table["COREX"];  
693 - const GFitsTableCol* ptr_corey = table["COREY"];  
694 - const GFitsTableCol* ptr_core_err = table["CORE_ERR"];  
695 - const GFitsTableCol* ptr_xmax = table["XMAX"];  
696 - const GFitsTableCol* ptr_xmax_err = table["XMAX_ERR"];  
697 - const GFitsTableCol* ptr_energy = table["ENERGY"];  
698 - const GFitsTableCol* ptr_energy_err = table["ENERGY_ERR"];  
699 -  
700 - // Check for phase column  
701 - const GFitsTableCol* ptr_phase;  
702 - if (table.contains("PHASE")) {  
703 - m_has_phase = true;  
704 - ptr_phase = table["PHASE"]; 671 + const GFitsTableCol* ptr_eid = table["EVENT_ID"];
  672 + const GFitsTableCol* ptr_time = table["TIME"];
  673 + const GFitsTableCol* ptr_ra = table["RA"];
  674 + const GFitsTableCol* ptr_dec = table["DEC"];
  675 + const GFitsTableCol* ptr_energy = table["ENERGY"];
  676 +
  677 + // Check for DETX and DETY columns
  678 + const GFitsTableCol* ptr_detx;
  679 + const GFitsTableCol* ptr_dety;
  680 + if (table.contains("DETX") && table.contains("DETY")) {
  681 + m_has_detxy = true;
  682 + ptr_detx = table["DETX"];
  683 + ptr_dety = table["DETY"];
705 } 684 }
706 else { 685 else {
707 - m_has_phase = false;  
708 - }  
709 -  
710 - // Copy data from columns into GCTAEventAtom objects  
711 - GCTAEventAtom event;  
712 - for (int i = 0; i < num; ++i) {  
713 - event.m_index = i;  
714 - event.m_time.set(ptr_time->real(i), m_gti.reference());  
715 - event.m_dir.dir().radec_deg(ptr_ra->real(i), ptr_dec->real(i));  
716 - event.m_dir.detx(ptr_detx->real(i)*gammalib::deg2rad);  
717 - event.m_dir.dety(ptr_dety->real(i)*gammalib::deg2rad);  
718 - event.m_energy.TeV(ptr_energy->real(i));  
719 - event.m_event_id = ptr_eid->integer(i);  
720 - event.m_obs_id = 0;  
721 - event.m_multip = ptr_multip->integer(i);  
722 - event.m_telmask = 0;  
723 - event.m_dir_err = ptr_dir_err->real(i);  
724 - event.m_alt = ptr_alt->real(i);  
725 - event.m_az = ptr_az->real(i);  
726 - event.m_corex = ptr_corex->real(i);  
727 - event.m_corey = ptr_corey->real(i);  
728 - event.m_core_err = ptr_core_err->real(i);  
729 - event.m_xmax = ptr_xmax->real(i);  
730 - event.m_xmax_err = ptr_xmax_err->real(i);  
731 - event.m_shwidth = 0.0;  
732 - event.m_shlength = 0.0;  
733 - event.m_energy_err = ptr_energy_err->real(i);  
734 -  
735 - // Set pulse phase if available  
736 - if (m_has_phase) {  
737 - event.m_phase = ptr_phase->real(i);  
738 - }  
739 -  
740 - m_events.push_back(event); 686 + m_has_detxy = false;
741 } 687 }
742 688
743 - } // endif: there were events  
744 -  
745 - // Return  
746 - return;  
747 -}  
748 -  
749 -  
750 -/***********************************************************************//**  
751 - * @brief Read CTA events from FITS table (version 1)  
752 - *  
753 - * @param[in] table FITS table.  
754 - *  
755 - * This method reads the CTA event list from a FITS table HDU into memory.  
756 - *  
757 - * @todo Implement agreed column format  
758 - ***************************************************************************/  
759 -void GCTAEventList::read_events_v1(const GFitsTable& table)  
760 -{  
761 - // Clear existing events  
762 - m_events.clear();  
763 -  
764 - // Extract number of events in FITS file  
765 - int num = table.nrows();  
766 -  
767 - // If there are events then load them  
768 - if (num > 0) {  
769 -  
770 - // Reserve data  
771 - m_events.reserve(num);  
772 -  
773 - // Get column pointers  
774 - const GFitsTableCol* ptr_eid = table["EVENT_ID"];  
775 - const GFitsTableCol* ptr_oid = table["OBS_ID"];  
776 - const GFitsTableCol* ptr_time = table["TIME"];  
777 - const GFitsTableCol* ptr_multip = table["MULTIP"];  
778 - const GFitsTableCol* ptr_ra = table["RA"];  
779 - const GFitsTableCol* ptr_dec = table["DEC"];  
780 - const GFitsTableCol* ptr_dir_err = table["DIR_ERR"];  
781 - const GFitsTableCol* ptr_detx = table["DETX"];  
782 - const GFitsTableCol* ptr_dety = table["DETY"];  
783 - const GFitsTableCol* ptr_alt = table["ALT"];  
784 - const GFitsTableCol* ptr_az = table["AZ"];  
785 - const GFitsTableCol* ptr_corex = table["COREX"];  
786 - const GFitsTableCol* ptr_corey = table["COREY"];  
787 - const GFitsTableCol* ptr_core_err = table["CORE_ERR"];  
788 - const GFitsTableCol* ptr_xmax = table["XMAX"];  
789 - const GFitsTableCol* ptr_xmax_err = table["XMAX_ERR"];  
790 - const GFitsTableCol* ptr_shw = table["SHWIDTH"];  
791 - const GFitsTableCol* ptr_shl = table["SHLENGTH"];  
792 - const GFitsTableCol* ptr_energy = table["ENERGY"];  
793 - const GFitsTableCol* ptr_energy_err = table["ENERGY_ERR"];  
794 -  
795 - // Check for phase column 689 + // Check for PHASE column
796 const GFitsTableCol* ptr_phase; 690 const GFitsTableCol* ptr_phase;
797 if (table.contains("PHASE")) { 691 if (table.contains("PHASE")) {
798 m_has_phase = true; 692 m_has_phase = true;
@@ -805,98 +699,57 @@ void GCTAEventList::read_events_v1(const GFitsTable&amp; table) @@ -805,98 +699,57 @@ void GCTAEventList::read_events_v1(const GFitsTable&amp; table)
805 // Copy data from columns into GCTAEventAtom objects 699 // Copy data from columns into GCTAEventAtom objects
806 GCTAEventAtom event; 700 GCTAEventAtom event;
807 for (int i = 0; i < num; ++i) { 701 for (int i = 0; i < num; ++i) {
808 - event.m_index = i; 702 +
  703 + // Set mandatory information
809 event.m_time.set(ptr_time->real(i), m_gti.reference()); 704 event.m_time.set(ptr_time->real(i), m_gti.reference());
810 event.m_dir.dir().radec_deg(ptr_ra->real(i), ptr_dec->real(i)); 705 event.m_dir.dir().radec_deg(ptr_ra->real(i), ptr_dec->real(i));
811 - event.m_dir.detx(ptr_detx->real(i)*gammalib::deg2rad);  
812 - event.m_dir.dety(ptr_dety->real(i)*gammalib::deg2rad);  
813 event.m_energy.TeV(ptr_energy->real(i)); 706 event.m_energy.TeV(ptr_energy->real(i));
814 - event.m_event_id = ptr_eid->integer(i);  
815 - event.m_obs_id = ptr_oid->integer(i);  
816 - event.m_multip = ptr_multip->integer(i);  
817 - event.m_telmask = 0;  
818 - event.m_dir_err = ptr_dir_err->real(i);  
819 - event.m_alt = ptr_alt->real(i);  
820 - event.m_az = ptr_az->real(i);  
821 - event.m_corex = ptr_corex->real(i);  
822 - event.m_corey = ptr_corey->real(i);  
823 - event.m_core_err = ptr_core_err->real(i);  
824 - event.m_xmax = ptr_xmax->real(i);  
825 - event.m_xmax_err = ptr_xmax_err->real(i);  
826 - event.m_shwidth = ptr_shw->real(i);  
827 - event.m_shlength = ptr_shl->real(i);  
828 - event.m_energy_err = ptr_energy_err->real(i);  
829 -  
830 - // Set phase if available 707 + event.m_event_id = ptr_eid->integer(i);
  708 + event.m_index = i;
  709 +
  710 + // If available, set detector coordinates in radians
  711 + if (m_has_detxy) {
  712 + event.m_dir.detx(ptr_detx->real(i)*gammalib::deg2rad);
  713 + event.m_dir.dety(ptr_dety->real(i)*gammalib::deg2rad);
  714 + }
  715 +
  716 + // If available, set pulse phase
831 if (m_has_phase) { 717 if (m_has_phase) {
832 event.m_phase = ptr_phase->real(i); 718 event.m_phase = ptr_phase->real(i);
833 } 719 }
834 720
  721 + // Append event
835 m_events.push_back(event); 722 m_events.push_back(event);
836 - }  
837 723
838 - } // endif: there were events 724 + } // endfor: looped over all events
839 725
840 - // Return  
841 - return;  
842 -} 726 + // Loop over table and find optional columns
  727 + for (int i = 0; i < table.ncols(); ++i) {
843 728
  729 + // Get column pointer
  730 + const GFitsTableCol* col = table[i];
844 731
845 -/***********************************************************************//**  
846 - * @brief Read Hillas information for CTA events from FITS table  
847 - *  
848 - * @param[in] table FITS table.  
849 - *  
850 - * This method reads the Hillas reconstruction information for CTA events  
851 - * from an EVENTS file. It searches for the columns HIL_MSW, HIL_MSW_ERR,  
852 - * HIL_MSL, and HIL_MSL_ERR in the FITS table and extracts the relevant  
853 - * columns from the FITS file. If a column is not found, no action is  
854 - * performed.  
855 - *  
856 - * @todo Verify consistency of event list size  
857 - ***************************************************************************/  
858 -void GCTAEventList::read_events_hillas(const GFitsTable& table)  
859 -{  
860 - // Extract number of events in FITS file  
861 - int num = table.nrows(); 732 + // Get column name
  733 + const std::string name = col->name();
862 734
863 - //TODO: Make sure that dimension is consistent with existing event  
864 - // list  
865 -  
866 - // Continue only if there are events  
867 - if (num > 0) {  
868 -  
869 - // HIL_MSW  
870 - if (table.contains("HIL_MSW")) {  
871 - const GFitsTableCol* ptr = table["HIL_MSW"];  
872 - for (int i = 0; i < num; ++i) {  
873 - m_events[i].m_hil_msw = ptr->real(i);  
874 - }  
875 - }  
876 -  
877 - // HIL_MSW_ERR  
878 - if (table.contains("HIL_MSW_ERR")) {  
879 - const GFitsTableCol* ptr = table["HIL_MSW_ERR"];  
880 - for (int i = 0; i < num; ++i) {  
881 - m_events[i].m_hil_msw_err = ptr->real(i); 735 + // If column was mandatory or optional then skip it ...
  736 + if (name == "EVENT_ID" ||
  737 + name == "TIME" ||
  738 + name == "RA" ||
  739 + name == "DEC" ||
  740 + name == "ENERGY" ||
  741 + name == "DETX" ||
  742 + name == "DETY" ||
  743 + name == "PHASE") {
  744 + continue;
882 } 745 }
883 - }  
884 746
885 - // HIL_MSL  
886 - if (table.contains("HIL_MSL")) {  
887 - const GFitsTableCol* ptr = table["HIL_MSL"];  
888 - for (int i = 0; i < num; ++i) {  
889 - m_events[i].m_hil_msl = ptr->real(i); 747 + // ... otherwise keep a copy
  748 + else {
  749 + m_columns.push_back(col->clone());
890 } 750 }
891 - }  
892 751
893 - // HIL_MSL_ERR  
894 - if (table.contains("HIL_MSL_ERR")) {  
895 - const GFitsTableCol* ptr = table["HIL_MSL_ERR"];  
896 - for (int i = 0; i < num; ++i) {  
897 - m_events[i].m_hil_msl_err = ptr->real(i);  
898 - }  
899 - } 752 + } // endfor: loop over optional columns
900 753
901 } // endif: there were events 754 } // endif: there were events
902 755
@@ -908,12 +761,26 @@ void GCTAEventList::read_events_hillas(const GFitsTable&amp; table) @@ -908,12 +761,26 @@ void GCTAEventList::read_events_hillas(const GFitsTable&amp; table)
908 /***********************************************************************//** 761 /***********************************************************************//**
909 * @brief Write CTA events into FITS table 762 * @brief Write CTA events into FITS table
910 * 763 *
911 - * @param[in] hdu FITS table HDU. 764 + * @param[in] table FITS binary table.
912 * 765 *
913 - * Write the CTA event list into FITS table. 766 + * Write the CTA event list into a FITS binary table.
914 * 767 *
915 - * @todo The TELMASK column is allocated with a dummy length of 100.  
916 - * @todo Implement agreed column format 768 + * The following mandatory columns will be written into the FITS table:
  769 + *
  770 + * - EVENT_ID
  771 + * - TIME
  772 + * - RA
  773 + * - DEC
  774 + * - ENERGY
  775 + *
  776 + * If available, also the
  777 + *
  778 + * - DETX
  779 + * - DETY
  780 + * - PHASE
  781 + *
  782 + * columns are written. Any other columns that were read by the read_events()
  783 + * method will be also written into the table.
917 ***************************************************************************/ 784 ***************************************************************************/
918 void GCTAEventList::write_events(GFitsBinTable& hdu) const 785 void GCTAEventList::write_events(GFitsBinTable& hdu) const
919 { 786 {
@@ -923,109 +790,55 @@ void GCTAEventList::write_events(GFitsBinTable&amp; hdu) const @@ -923,109 +790,55 @@ void GCTAEventList::write_events(GFitsBinTable&amp; hdu) const
923 // If there are events then write them now 790 // If there are events then write them now
924 if (size() > 0) { 791 if (size() > 0) {
925 792
926 - // Allocate columns  
927 - GFitsTableULongCol col_eid = GFitsTableULongCol("EVENT_ID", size());  
928 - GFitsTableULongCol col_oid = GFitsTableULongCol("OBS_ID", size());  
929 - GFitsTableDoubleCol col_time = GFitsTableDoubleCol("TIME", size());  
930 - GFitsTableDoubleCol col_live = GFitsTableDoubleCol("TLIVE", size());  
931 - GFitsTableShortCol col_multip = GFitsTableShortCol("MULTIP", size());  
932 - GFitsTableBitCol col_telmask = GFitsTableBitCol("TELMASK", size(), 100);  
933 - GFitsTableFloatCol col_ra = GFitsTableFloatCol("RA", size());  
934 - GFitsTableFloatCol col_dec = GFitsTableFloatCol("DEC", size());  
935 - GFitsTableFloatCol col_direrr = GFitsTableFloatCol("DIR_ERR", size());  
936 - GFitsTableFloatCol col_detx = GFitsTableFloatCol("DETX", size());  
937 - GFitsTableFloatCol col_dety = GFitsTableFloatCol("DETY", size());  
938 - GFitsTableFloatCol col_alt = GFitsTableFloatCol("ALT", size());  
939 - GFitsTableFloatCol col_az = GFitsTableFloatCol("AZ", size());  
940 - GFitsTableFloatCol col_corex = GFitsTableFloatCol("COREX", size());  
941 - GFitsTableFloatCol col_corey = GFitsTableFloatCol("COREY", size());  
942 - GFitsTableFloatCol col_core_err = GFitsTableFloatCol("CORE_ERR", size());  
943 - GFitsTableFloatCol col_xmax = GFitsTableFloatCol("XMAX", size());  
944 - GFitsTableFloatCol col_xmax_err = GFitsTableFloatCol("XMAX_ERR", size());  
945 - GFitsTableFloatCol col_shw = GFitsTableFloatCol("SHWIDTH", size());  
946 - GFitsTableFloatCol col_shl = GFitsTableFloatCol("SHLENGTH", size());  
947 - GFitsTableFloatCol col_energy = GFitsTableFloatCol("ENERGY", size());  
948 - GFitsTableFloatCol col_energy_err = GFitsTableFloatCol("ENERGY_ERR", size());  
949 - GFitsTableFloatCol col_hil_msw = GFitsTableFloatCol("HIL_MSW", size());  
950 - GFitsTableFloatCol col_hil_msw_err = GFitsTableFloatCol("HIL_MSW_ERR", size());  
951 - GFitsTableFloatCol col_hil_msl = GFitsTableFloatCol("HIL_MSL", size());  
952 - GFitsTableFloatCol col_hil_msl_err = GFitsTableFloatCol("HIL_MSL_ERR", size());  
953 -  
954 - // Create phase column if required  
955 - GFitsTableFloatCol col_phase;  
956 - if (m_has_phase) {  
957 - col_phase = GFitsTableFloatCol("PHASE", size());  
958 - } 793 + // Allocate mandatory columns
  794 + GFitsTableULongCol col_eid = GFitsTableULongCol("EVENT_ID", size());
  795 + GFitsTableDoubleCol col_time = GFitsTableDoubleCol("TIME", size());
  796 + GFitsTableFloatCol col_ra = GFitsTableFloatCol("RA", size());
  797 + GFitsTableFloatCol col_dec = GFitsTableFloatCol("DEC", size());
  798 + GFitsTableFloatCol col_energy = GFitsTableFloatCol("ENERGY", size());
959 799
960 - // Fill columns 800 + // Fill mandatory columns
961 for (int i = 0; i < size(); ++i) { 801 for (int i = 0; i < size(); ++i) {
962 - col_eid(i) = m_events[i].m_event_id;  
963 - col_oid(i) = m_events[i].m_obs_id;  
964 - col_time(i) = m_events[i].time().convert(m_gti.reference());  
965 - col_live(i) = 0.0;  
966 - col_multip(i) = 0;  
967 - //col_telmask  
968 - col_ra(i) = m_events[i].dir().dir().ra_deg();  
969 - col_dec(i) = m_events[i].dir().dir().dec_deg();  
970 - col_direrr(i) = m_events[i].m_dir_err;  
971 - col_detx(i) = m_events[i].dir().detx() * gammalib::rad2deg;  
972 - col_dety(i) = m_events[i].dir().dety() * gammalib::rad2deg;  
973 - col_alt(i) = m_events[i].m_alt;  
974 - col_az(i) = m_events[i].m_az;  
975 - col_corex(i) = m_events[i].m_corex;  
976 - col_corey(i) = m_events[i].m_corey;  
977 - col_core_err(i) = m_events[i].m_core_err;  
978 - col_xmax(i) = m_events[i].m_xmax;  
979 - col_xmax_err(i) = m_events[i].m_xmax_err;  
980 - col_shw(i) = m_events[i].m_shwidth;  
981 - col_shl(i) = m_events[i].m_shlength;  
982 - col_energy(i) = m_events[i].energy().TeV();  
983 - col_energy_err(i) = m_events[i].m_energy_err;  
984 - col_hil_msw(i) = m_events[i].m_hil_msw;  
985 - col_hil_msw_err(i) = m_events[i].m_hil_msw_err;  
986 - col_hil_msl(i) = m_events[i].m_hil_msl;  
987 - col_hil_msl_err(i) = m_events[i].m_hil_msl_err;  
988 -  
989 - // Optionally fill pulse phase column  
990 - if (m_has_phase) {  
991 - col_phase(i) = m_events[i].m_phase;  
992 - }  
993 -  
994 - } // endfor: looped over rows 802 + col_eid(i) = m_events[i].m_event_id;
  803 + col_time(i) = m_events[i].time().convert(m_gti.reference());
  804 + col_ra(i) = m_events[i].dir().dir().ra_deg();
  805 + col_dec(i) = m_events[i].dir().dir().dec_deg();
  806 + col_energy(i) = m_events[i].energy().TeV();
  807 + }
995 808
996 - // Append columns to table 809 + // Append mandatory columns to table
997 hdu.append(col_eid); 810 hdu.append(col_eid);
998 - hdu.append(col_oid);  
999 hdu.append(col_time); 811 hdu.append(col_time);
1000 - hdu.append(col_live);  
1001 - hdu.append(col_multip);  
1002 - hdu.append(col_telmask);  
1003 hdu.append(col_ra); 812 hdu.append(col_ra);
1004 hdu.append(col_dec); 813 hdu.append(col_dec);
1005 - hdu.append(col_direrr);  
1006 - hdu.append(col_detx);  
1007 - hdu.append(col_dety);  
1008 - hdu.append(col_alt);  
1009 - hdu.append(col_az);  
1010 - hdu.append(col_corex);  
1011 - hdu.append(col_corey);  
1012 - hdu.append(col_core_err);  
1013 - hdu.append(col_xmax);  
1014 - hdu.append(col_xmax_err);  
1015 - hdu.append(col_shw);  
1016 - hdu.append(col_shl);  
1017 hdu.append(col_energy); 814 hdu.append(col_energy);
1018 - hdu.append(col_energy_err);  
1019 - hdu.append(col_hil_msw);  
1020 - hdu.append(col_hil_msw_err);  
1021 - hdu.append(col_hil_msl);  
1022 - hdu.append(col_hil_msl_err);  
1023 815
1024 - // Optionally write phase column 816 + // If available, add detector coordinates in degrees
  817 + if (m_has_detxy) {
  818 + GFitsTableFloatCol col_detx = GFitsTableFloatCol("DETX", size());
  819 + GFitsTableFloatCol col_dety = GFitsTableFloatCol("DETY", size());
  820 + for (int i = 0; i < size(); ++i) {
  821 + col_detx(i) = m_events[i].dir().detx() * gammalib::rad2deg;
  822 + col_dety(i) = m_events[i].dir().dety() * gammalib::rad2deg;
  823 + }
  824 + hdu.append(col_detx);
  825 + hdu.append(col_dety);
  826 + }
  827 +
  828 + // If available, add event phase
1025 if (m_has_phase) { 829 if (m_has_phase) {
  830 + GFitsTableFloatCol col_phase = GFitsTableFloatCol("PHASE", size());
  831 + for (int i = 0; i < size(); ++i) {
  832 + col_phase(i) = m_events[i].m_phase;
  833 + }
1026 hdu.append(col_phase); 834 hdu.append(col_phase);
1027 } 835 }
1028 836
  837 + // Append other columns to table
  838 + for (int i = 0; i < m_columns.size(); ++i) {
  839 + hdu.append(*m_columns[i]);
  840 + }
  841 +
1029 } // endif: there were events to write 842 } // endif: there were events to write
1030 843
1031 // Return 844 // Return