Domenico Tiziani / gammalib

Commit 053297644ceead80d691973e662c27ff86e9c95c

Authored by Michael Mayer
1 parent 1372c6c5

Allow reading/writing of necessary and optional columns for GCTAEventList and GCTAEventAtom

Inline Side-by-side
Showing 5 changed files with 82 additions and 385 deletions
inst/cta/include/GCTAEventAtom.hpp
  View file @0532976
@@ -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 @0532976
@@ -89,9 +89,6 @@ protected: @@ -89,9 +89,6 @@ protected:
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); 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; 92 void write_events(GFitsBinTable& hdu) 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;
@@ -100,11 +97,15 @@ protected: @@ -100,11 +97,15 @@ protected:
100 // Protected members 97 // Protected members
101 GCTARoi m_roi; //!< Region of interest 98 GCTARoi m_roi; //!< Region of interest
102 std::vector<GCTAEventAtom> m_events; //!< Events 99 std::vector<GCTAEventAtom> m_events; //!< Events
  100 + std::vector<GFitsTableCol*> m_columns; //!< Pointers to optional columns
103 bool m_has_phase; //!< Signal presence of phase 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
107 mutable std::vector<std::vector<double> > m_irf_values; //!< IRF values 106 mutable std::vector<std::vector<double> > m_irf_values; //!< IRF values
  107 +
  108 +
108 }; 109 };
109 110
110 111
inst/cta/pyext/GCTAEventAtom.i
  View file @0532976
@@ -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 @0532976
@@ -209,24 +209,6 @@ void GCTAEventAtom::init_members(void) @@ -209,24 +209,6 @@ void GCTAEventAtom::init_members(void)
209 m_energy.clear(); 209 m_energy.clear();
210 m_index = 0; 210 m_index = 0;
211 m_event_id = 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; 212 m_phase = 0.0;
231 213
232 // Return 214 // Return
@@ -247,24 +229,6 @@ void GCTAEventAtom::copy_members(const GCTAEventAtom&amp; atom) @@ -247,24 +229,6 @@ void GCTAEventAtom::copy_members(const GCTAEventAtom&amp; atom)
247 m_energy = atom.m_energy; 229 m_energy = atom.m_energy;
248 m_index = atom.m_index; 230 m_index = atom.m_index;
249 m_event_id = atom.m_event_id; 231 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; 232 m_phase = atom.m_phase;
269 233
270 // Return 234 // Return
inst/cta/src/GCTAEventList.cpp
  View file @0532976
@@ -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 }
@@ -620,8 +635,9 @@ void GCTAEventList::free_members(void) @@ -620,8 +635,9 @@ void GCTAEventList::free_members(void)
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. 637 * 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. 638 + * It reads required columns for analysis such as "EVENT_ID", "OBS_ID",
  639 + * "RA", "DEC", "TIME" and "ENERGY". All other columns are stored as well
  640 + * and will be written out in the write()-method.
625 ***************************************************************************/ 641 ***************************************************************************/
626 void GCTAEventList::read_events(const GFitsTable& table) 642 void GCTAEventList::read_events(const GFitsTable& table)
627 { 643 {
@@ -629,47 +645,6 @@ void GCTAEventList::read_events(const GFitsTable&amp; table) @@ -629,47 +645,6 @@ void GCTAEventList::read_events(const GFitsTable&amp; table)
629 m_events.clear(); 645 m_events.clear();
630 646
631 // Extract number of events in FITS file 647 // 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)  
659 - *  
660 - * @param[in] table FITS table.  
661 - *  
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.  
666 - ***************************************************************************/  
667 -void GCTAEventList::read_events_v0(const GFitsTable& table)  
668 -{  
669 - // Clear existing events  
670 - m_events.clear();  
671 -  
672 - // Extract number of events in FITS file  
673 int num = table.nrows(); 648 int num = table.nrows();
674 649
675 // If there are events then load them 650 // If there are events then load them
@@ -681,21 +656,9 @@ void GCTAEventList::read_events_v0(const GFitsTable&amp; table) @@ -681,21 +656,9 @@ void GCTAEventList::read_events_v0(const GFitsTable&amp; table)
681 // Get column pointers 656 // Get column pointers
682 const GFitsTableCol* ptr_eid = table["EVENT_ID"]; 657 const GFitsTableCol* ptr_eid = table["EVENT_ID"];
683 const GFitsTableCol* ptr_time = table["TIME"]; 658 const GFitsTableCol* ptr_time = table["TIME"];
684 - const GFitsTableCol* ptr_multip = table["MULTIP"];  
685 const GFitsTableCol* ptr_ra = table["RA"]; 659 const GFitsTableCol* ptr_ra = table["RA"];
686 const GFitsTableCol* ptr_dec = table["DEC"]; 660 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"]; 661 const GFitsTableCol* ptr_energy = table["ENERGY"];
698 - const GFitsTableCol* ptr_energy_err = table["ENERGY_ERR"];  
699 662
700 // Check for phase column 663 // Check for phase column
701 const GFitsTableCol* ptr_phase; 664 const GFitsTableCol* ptr_phase;
@@ -707,99 +670,16 @@ void GCTAEventList::read_events_v0(const GFitsTable&amp; table) @@ -707,99 +670,16 @@ void GCTAEventList::read_events_v0(const GFitsTable&amp; table)
707 m_has_phase = false; 670 m_has_phase = false;
708 } 671 }
709 672
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);  
741 - }  
742 -  
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  
796 - const GFitsTableCol* ptr_phase;  
797 - if (table.contains("PHASE")) {  
798 - m_has_phase = true;  
799 - ptr_phase = table["PHASE"]; 673 + // Check for DETX and DETY columns
  674 + const GFitsTableCol* ptr_detx;
  675 + const GFitsTableCol* ptr_dety;
  676 + if (table.contains("DETX") && table.contains("DETY")) {
  677 + m_has_detxy = true;
  678 + ptr_detx = table["DETX"];
  679 + ptr_dety = table["DETY"];
800 } 680 }
801 else { 681 else {
802 - m_has_phase = false; 682 + m_has_detxy = false;
803 } 683 }
804 684
805 // Copy data from columns into GCTAEventAtom objects 685 // Copy data from columns into GCTAEventAtom objects
@@ -808,95 +688,45 @@ void GCTAEventList::read_events_v1(const GFitsTable&amp; table) @@ -808,95 +688,45 @@ void GCTAEventList::read_events_v1(const GFitsTable&amp; table)
808 event.m_index = i; 688 event.m_index = i;
809 event.m_time.set(ptr_time->real(i), m_gti.reference()); 689 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)); 690 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)); 691 event.m_energy.TeV(ptr_energy->real(i));
814 event.m_event_id = ptr_eid->integer(i); 692 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 693 +
  694 + // Set detector coordinates if available
  695 + if (m_has_detxy) {
  696 + event.m_dir.detx(ptr_detx->real(i)*gammalib::deg2rad);
  697 + event.m_dir.dety(ptr_dety->real(i)*gammalib::deg2rad);
  698 + }
  699 +
  700 + // Set pulse phase if available
831 if (m_has_phase) { 701 if (m_has_phase) {
832 event.m_phase = ptr_phase->real(i); 702 event.m_phase = ptr_phase->real(i);
833 } 703 }
834 704
  705 + // Append event
835 m_events.push_back(event); 706 m_events.push_back(event);
836 } 707 }
837 708
838 - } // endif: there were events  
839 -  
840 - // Return  
841 - return;  
842 -}  
843 -  
844 -  
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();  
862 -  
863 - //TODO: Make sure that dimension is consistent with existing event  
864 - // list 709 + // Loop over table and find optional columns
  710 + for (int i = 0; i < table.ncols(); ++i) {
865 711
866 - // Continue only if there are events  
867 - if (num > 0) { 712 + // Get column pointer
  713 + const GFitsTableCol* col = table[i];
868 714
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 - } 715 + // Get column name
  716 + const std::string name = col->name();
876 717
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); 718 + // Check if column was handled before
  719 + if (name == "EVENT_ID" || name == "ENERGY" ||
  720 + name == "RA" || name == "DEC" ||
  721 + name == "TIME" || name == "DETX" ||
  722 + name == "DETY" || name == "PHASE") {
  723 + continue;
882 } 724 }
883 - }  
884 -  
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); 725 + else {
  726 + m_columns.push_back(col->clone());
890 } 727 }
891 - }  
892 728
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 - } 729 + } // endfor: loop over optional columns
900 730
901 } // endif: there were events 731 } // endif: there were events
902 732
@@ -912,7 +742,6 @@ void GCTAEventList::read_events_hillas(const GFitsTable&amp; table) @@ -912,7 +742,6 @@ void GCTAEventList::read_events_hillas(const GFitsTable&amp; table)
912 * 742 *
913 * Write the CTA event list into FITS table. 743 * Write the CTA event list into FITS table.
914 * 744 *
915 - * @todo The TELMASK column is allocated with a dummy length of 100.  
916 * @todo Implement agreed column format 745 * @todo Implement agreed column format
917 ***************************************************************************/ 746 ***************************************************************************/
918 void GCTAEventList::write_events(GFitsBinTable& hdu) const 747 void GCTAEventList::write_events(GFitsBinTable& hdu) const
@@ -923,109 +752,59 @@ void GCTAEventList::write_events(GFitsBinTable&amp; hdu) const @@ -923,109 +752,59 @@ void GCTAEventList::write_events(GFitsBinTable&amp; hdu) const
923 // If there are events then write them now 752 // If there are events then write them now
924 if (size() > 0) { 753 if (size() > 0) {
925 754
926 - // Allocate columns 755 + // Allocate required columns
927 GFitsTableULongCol col_eid = GFitsTableULongCol("EVENT_ID", size()); 756 GFitsTableULongCol col_eid = GFitsTableULongCol("EVENT_ID", size());
928 - GFitsTableULongCol col_oid = GFitsTableULongCol("OBS_ID", size());  
929 GFitsTableDoubleCol col_time = GFitsTableDoubleCol("TIME", size()); 757 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()); 758 GFitsTableFloatCol col_ra = GFitsTableFloatCol("RA", size());
934 GFitsTableFloatCol col_dec = GFitsTableFloatCol("DEC", size()); 759 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()); 760 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 - }  
959 761
960 // Fill columns 762 // Fill columns
961 for (int i = 0; i < size(); ++i) { 763 for (int i = 0; i < size(); ++i) {
962 col_eid(i) = m_events[i].m_event_id; 764 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()); 765 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(); 766 col_ra(i) = m_events[i].dir().dir().ra_deg();
969 col_dec(i) = m_events[i].dir().dir().dec_deg(); 767 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(); 768 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 769 + } // endfor: loop over columns
995 770
996 // Append columns to table 771 // Append columns to table
997 hdu.append(col_eid); 772 hdu.append(col_eid);
998 - hdu.append(col_oid);  
999 hdu.append(col_time); 773 hdu.append(col_time);
1000 - hdu.append(col_live);  
1001 - hdu.append(col_multip);  
1002 - hdu.append(col_telmask);  
1003 hdu.append(col_ra); 774 hdu.append(col_ra);
1004 hdu.append(col_dec); 775 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); 776 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 777
1024 - // Optionally write phase column 778 + // Fill detector coordinates if required
  779 + if (m_has_detxy) {
  780 + GFitsTableFloatCol col_detx = GFitsTableFloatCol("DETX", size());
  781 + GFitsTableFloatCol col_dety = GFitsTableFloatCol("DETY", size());
  782 + for (int i = 0; i < size(); ++i) {
  783 + col_detx(i) = m_events[i].dir().detx() * gammalib::rad2deg;
  784 + col_dety(i) = m_events[i].dir().dety() * gammalib::rad2deg;
  785 + }
  786 +
  787 + // Append columns to table
  788 + hdu.append(col_detx);
  789 + hdu.append(col_dety);
  790 + }
  791 +
  792 + // Fill pulse phase if required
1025 if (m_has_phase) { 793 if (m_has_phase) {
  794 + GFitsTableFloatCol col_phase = GFitsTableFloatCol("PHASE", size());
  795 + for (int i = 0; i < size(); ++i) {
  796 + col_phase(i) = m_events[i].m_phase;
  797 + }
  798 +
  799 + // Append column to table
1026 hdu.append(col_phase); 800 hdu.append(col_phase);
1027 } 801 }
1028 802
  803 + // Append other columns to table
  804 + for (int i = 0; i < m_columns.size(); ++i) {
  805 + hdu.append(*m_columns[i]);
  806 + }
  807 +
1029 } // endif: there were events to write 808 } // endif: there were events to write
1030 809
1031 // Return 810 // Return