GammaLib  2.0.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GEbounds.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GEbounds.cpp - Energy boundary class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2009-2020 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GEbounds.cpp
23  * @brief Energy boundary class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GFilename.hpp"
35 #include "GEbounds.hpp"
36 #include "GEnergies.hpp"
37 #include "GFits.hpp"
38 #include "GFitsTable.hpp"
39 #include "GFitsBinTable.hpp"
40 #include "GFitsTableDoubleCol.hpp"
41 #include "GXmlElement.hpp"
42 
43 /* __ Method name definitions ____________________________________________ */
44 #define G_READ_XML "GEbounds::read(GXmlElement&)"
45 #define G_WRITE_XML "GEbounds::write(GXmlElement&)"
46 #define G_REMOVE "GEbounds::remove(int&)"
47 #define G_EMIN_SET "GEbounds::emin(int&, GEnergy&)"
48 #define G_EMAX_SET "GEbounds::emax(int&, GEnergy&)"
49 #define G_EMIN_GET "GEbounds::emin(int&)"
50 #define G_EMAX_GET "GEbounds::emax(int&)"
51 #define G_EMEAN "GEbounds::emean(int&)"
52 #define G_ELOGMEAN "GEbounds::elogmean(int&)"
53 #define G_EWIDTH "GEbounds::ewidth(int&)"
54 #define G_INSERT_ENG "GEbounds::insert_eng(int&, GEnergy&, GEnergy&)"
55 
56 /* __ Macros _____________________________________________________________ */
57 
58 /* __ Coding definitions _________________________________________________ */
59 
60 /* __ Debug definitions __________________________________________________ */
61 
62 
63 /*==========================================================================
64  = =
65  = Constructors/destructors =
66  = =
67  ==========================================================================*/
68 
69 /***********************************************************************//**
70  * @brief Constructor
71  ***************************************************************************/
73 {
74  // Initialise class members for clean destruction
75  init_members();
76 
77  // Return
78  return;
79 }
80 
81 
82 /***********************************************************************//**
83  * @brief FITS file constructor
84  *
85  * @param[in] filename FITS file name.
86  *
87  * Constructs energy boundaries from a FITS file.
88  ***************************************************************************/
89 GEbounds::GEbounds(const GFilename& filename)
90 {
91  // Initialise members
92  init_members();
93 
94  // Load FITS file
95  load(filename);
96 
97  // Return
98  return;
99 }
100 
101 
102 /***********************************************************************//**
103  * @brief XML element constructor
104  *
105  * @param[in] xml XML element.
106  *
107  * Constructs energy boundaries from an XML element.
108  ***************************************************************************/
110 {
111  // Initialise members
112  init_members();
113 
114  // Read energy boundaries from XML element
115  read(xml);
116 
117  // Return
118  return;
119 }
120 
121 
122 /***********************************************************************//**
123  * @brief Energy container constructor
124  *
125  * @param[in] energies Energy container.
126  *
127  * Constructs energy boundaries from an energy container.
128  ***************************************************************************/
130 {
131  // Initialise members
132  init_members();
133 
134  // Set energy boundaries from energy container
135  set(energies);
136 
137  // Return
138  return;
139 }
140 
141 
142 /***********************************************************************//**
143  * @brief Copy constructor
144  *
145  * @param[in] ebds Energy boundaries.
146  ***************************************************************************/
148 {
149  // Initialise class members for clean destruction
150  init_members();
151 
152  // Copy members
153  copy_members(ebds);
154 
155  // Return
156  return;
157 }
158 
159 
160 /***********************************************************************//**
161  * @brief Single energy band constructor
162  *
163  * @param[in] emin Minimum energy of the interval.
164  * @param[in] emax Maximum energy of the interval.
165  *
166  * Constructs energy boundaries for one (emin, emax) energy band.
167  ***************************************************************************/
168 GEbounds::GEbounds(const GEnergy& emin, const GEnergy& emax)
169 {
170  // Initialise members
171  init_members();
172 
173  // Append energies
174  append(emin, emax);
175 
176  // Return
177  return;
178 }
179 
180 
181 /***********************************************************************//**
182  * @brief Interval constructor
183  *
184  * @param[in] num Number of energy intervals.
185  * @param[in] emin Minimum energy of first interval.
186  * @param[in] emax Maximum energy of last interval.
187  * @param[in] method Energy spacing method (one of "LIN", "LOG" or "POW").
188  * @param[in] gamma Power law index for @p POW method.
189  *
190  * Constructs energy boundaries by defining @p num successive energy
191  * intervals between @p emin and @p emax. The @p method parameter controls
192  * the energy spacing of the energy boundaries. See the set() method for
193  * more information.
194  ***************************************************************************/
195 GEbounds::GEbounds(const int& num,
196  const GEnergy& emin,
197  const GEnergy& emax,
198  const std::string& method,
199  const double& gamma)
200 {
201  // Initialise members
202  init_members();
203 
204  // Set intervals
205  set(num, emin, emax, method, gamma);
206 
207  // Return
208  return;
209 }
210 
211 
212 /***********************************************************************//**
213  * @brief Destructor
214  ***************************************************************************/
216 {
217  // Free members
218  free_members();
219 
220  // Return
221  return;
222 }
223 
224 
225 /*==========================================================================
226  = =
227  = Operators =
228  = =
229  ==========================================================================*/
230 
231 /***********************************************************************//**
232  * @brief Assignment operator
233  *
234  * @param[in] ebds Energy boundaries to be assigned.
235  * @return Energy boundaries.
236  ***************************************************************************/
238 {
239  // Execute only if object is not identical
240  if (this != &ebds) {
241 
242  // Free members
243  free_members();
244 
245  // Initialise private members for clean destruction
246  init_members();
247 
248  // Copy members
249  copy_members(ebds);
250 
251  } // endif: object was not identical
252 
253  // Return this object
254  return *this;
255 }
256 
257 
258 /*==========================================================================
259  = =
260  = Public methods =
261  = =
262  ==========================================================================*/
263 
264 /***********************************************************************//**
265  * @brief Clear energy boundaries
266  ***************************************************************************/
267 void GEbounds::clear(void)
268 {
269  // Free members
270  free_members();
271 
272  // Initialise private members
273  init_members();
274 
275  // Return
276  return;
277 }
278 
279 
280 /***********************************************************************//**
281  * @brief Clone energy boundaries
282  *
283  * @return Pointer to deep copy of energy boundaries.
284  ***************************************************************************/
286 {
287  return new GEbounds(*this);
288 }
289 
290 
291 /***********************************************************************//**
292  * @brief Append energy interval
293  *
294  * @param[in] emin Minimum energy of interval.
295  * @param[in] emax Maximum energy of interval.
296  *
297  * Appends an energy interval to the end of the energy boundaries
298  * container.
299  ***************************************************************************/
300 void GEbounds::append(const GEnergy& emin, const GEnergy& emax)
301 {
302  // Append interval at the end
303  insert_eng(m_num, emin, emax);
304 
305  // Return
306  return;
307 }
308 
309 
310 /***********************************************************************//**
311  * @brief Insert energy interval
312  *
313  * @param[in] emin Minimum energy of interval.
314  * @param[in] emax Maximum energy of interval.
315  *
316  * Inserts an energy interval into the energy boundaries after the first
317  * boundary that has a minimum energy smaller than @p emin. The method
318  * implicitely assumes that the intervals are ordered by increasing minimum
319  * energy.
320  ***************************************************************************/
321 void GEbounds::insert(const GEnergy& emin, const GEnergy& emax)
322 {
323  // Determine index at which interval should be inserted
324  int inx = 0;
325  for (; inx < m_num; ++inx) {
326  if (emin < m_min[inx])
327  break;
328  }
329 
330  // Insert interval
331  insert_eng(inx, emin, emax);
332 
333  // Return
334  return;
335 }
336 
337 
338 /***********************************************************************//**
339  * @brief Merge all overlapping or connecting successive energy intervals
340  *
341  * Merges all overlapping or connecting successive energy intervals. The
342  * method implicitely assumes that the intervals are ordered by increasing
343  * minimum energy.
344  *
345  * Note that the method does not actually reduce the memory size but just
346  * updates the information on the number of elements in the array.
347  ***************************************************************************/
348 void GEbounds::merge(void)
349 {
350  // Find overlaps
351  int i = 0;
352  int num = m_num;
353  while (i < num-1) {
354 
355  // If current energy interval overlaps with successor then merge both
356  // intervals, move all remaining intervals one position up, and
357  // reduce the number of elements
358  if (m_min[i+1] <= m_max[i]) {
359  m_min[i] = (m_min[i] < m_min[i+1]) ? m_min[i] : m_min[i+1];
360  m_max[i] = (m_max[i] > m_max[i+1]) ? m_max[i] : m_max[i+1];
361  for (int k = i+2; k < num; ++k) {
362  m_min[k-1] = m_min[k];
363  m_max[k-1] = m_max[k];
364  }
365  num--;
366  }
367 
368  // Otherwise increment interval index
369  else {
370  i++;
371  }
372 
373  } // endwhile: there were still intervals to check
374 
375  // Update number of elements in object
376  m_num = num;
377 
378  // Return
379  return;
380 }
381 
382 
383 /***********************************************************************//**
384  * @brief Merge energy interval into energy boundaries
385  *
386  * @param[in] emin Minimum energy of interval.
387  * @param[in] emax Maximum energy of interval.
388  *
389  * Inserts an energy interval into the energy boundaries after the first
390  * boundary that has a minimum energy smaller than @p emin and then merges
391  * any overlapping or connecting energy boundaries. The method implicitely
392  * assumes that the intervals are ordered by increasing minimum energy.
393  ***************************************************************************/
394 void GEbounds::merge(const GEnergy& emin, const GEnergy& emax)
395 {
396  // Determine index at which interval should be inserted
397  int inx = 0;
398  for (; inx < m_num; ++inx) {
399  if (emin < m_min[inx])
400  break;
401  }
402 
403  // Insert interval
404  insert_eng(inx, emin, emax);
405 
406  // Merge any overlapping energy intervals
407  merge();
408 
409  // Return
410  return;
411 }
412 
413 
414 /***********************************************************************//**
415  * @brief Remove energy interval
416  *
417  * @param[in] index Energy interval index (0,...,size()-1).
418  *
419  * Removes energy interval at @p index from the energy boundaries container.
420  * All intervals after the specified @p index are moved forward by one
421  * position.
422  *
423  * Note that the method does not actually reduce the memory size but just
424  * updates the information on the number of elements in the array.
425  ***************************************************************************/
426 void GEbounds::remove(const int& index)
427 {
428  #if defined(G_RANGE_CHECK)
429  // If index is outside boundary then throw an error
430  if (index < 0 || index >= m_num) {
431  throw GException::out_of_range(G_REMOVE, index, 0, m_num-1);
432  }
433  #endif
434 
435  // Move all elements located after index forward
436  for (int i = index+1; i < m_num; ++i) {
437  m_min[i-1] = m_min[i];
438  m_min[i-1] = m_max[i];
439  }
440 
441  // Reduce number of elements by one
442  m_num--;
443 
444  // Update attributes
445  set_attributes();
446 
447  // Return
448  return;
449 }
450 
451 
452 /***********************************************************************//**
453  * @brief Reserve space for energy intervals
454  *
455  * @param[in] num Number of elements.
456  *
457  * This method does nothing (it is needed to satisfy the GContainer
458  * interface).
459  ***************************************************************************/
460 void GEbounds::reserve(const int& num)
461 {
462  // Return
463  return;
464 }
465 
466 
467 /***********************************************************************//**
468  * @brief Append energy boundaries
469  *
470  * @param[in] ebds Energy boundaries.
471  *
472  * Append energy boundaries to the container.
473  ***************************************************************************/
474 void GEbounds::extend(const GEbounds& ebds)
475 {
476  // Do nothing if energy boundaries are empty
477  if (!ebds.is_empty()) {
478 
479  // Allocate new intervals
480  int num = m_num+ebds.size();
481  GEnergy* min = new GEnergy[num];
482  GEnergy* max = new GEnergy[num];
483 
484  // Initialise index
485  int inx = 0;
486 
487  // Copy existing intervals
488  for (; inx < m_num; ++inx) {
489  min[inx] = m_min[inx];
490  max[inx] = m_max[inx];
491  }
492 
493  // Append intervals
494  for (int i = 0; i < ebds.size(); ++i, ++inx) {
495  min[inx] = ebds.m_min[i];
496  max[inx] = ebds.m_max[i];
497  }
498 
499  // Free memory
500  if (m_min != NULL) delete [] m_min;
501  if (m_max != NULL) delete [] m_max;
502 
503  // Set new memory
504  m_min = min;
505  m_max = max;
506 
507  // Set number of elements
508  m_num = num;
509 
510  // Set attributes
511  set_attributes();
512 
513  } // endif: energy boundaries were not empty
514 
515  // Return
516  return;
517 }
518 
519 
520 /***********************************************************************//**
521  * @brief Set energy boundaries from energy container
522  *
523  * @param[in] energies Energy container.
524  *
525  * Sets the energy boundaries from an energy container. Each two subsequent
526  * energies in the energy container will form an energy boundary. This
527  * means that n energies will lead to n-1 energy boundaries with the
528  * following mapping:
529  *
530  * [energies[0], energies[1]]
531  * [energies[1], energies[2]]
532  * ...
533  * [energies[n-2], energies[n-1]]
534  *
535  * If there is only one energy in the container the following empty energy
536  * boundary will be appended:
537  *
538  * [energies[0], energies[0]]
539  ***************************************************************************/
540 void GEbounds::set(const GEnergies& energies)
541 {
542  // Initialise members
543  clear();
544 
545  // Get number of energies in container
546  int num = energies.size();
547 
548  // If there is only one energy in the container then append an empty
549  // energy boundary
550  if (num == 1) {
551  append(energies[0], energies[0]);
552  }
553 
554  // ... otherwise if there is more than one energy in the container
555  // then append subsequent energies as boundaries
556  else if (num > 1) {
557  for (int i = 0; i < num-1; ++i) {
558  append(energies[i], energies[i+1]);
559  }
560  }
561 
562  // Set attributes
563  set_attributes();
564 
565  // Return
566  return;
567 }
568 
569 
570 /***********************************************************************//**
571  * @brief Set energy intervals
572  *
573  * @param[in] num Number of energy intervals.
574  * @param[in] emin Minimum energy of first interval.
575  * @param[in] emax Maximum energy of last interval.
576  * @param[in] method Energy spacing method (one of "LIN", "LOG" or "POW").
577  * @param[in] gamma Power law index for @p POW method.
578  *
579  * Sets energy boundaries by defining @p num successive energy intervals
580  * between @p emin and @p emax. The @p method parameter controls the energy
581  * spacing of the energy boundaries. See the GEnergies::set() method for
582  * more information.
583  ***************************************************************************/
584 void GEbounds::set(const int& num,
585  const GEnergy& emin,
586  const GEnergy& emax,
587  const std::string& method,
588  const double& gamma)
589 {
590  // Set energies
591  GEnergies energies(num+1, emin, emax, method, gamma);
592 
593  // Set energy boundaries from energies
594  set(energies);
595 
596  // Return
597  return;
598 }
599 
600 
601 /***********************************************************************//**
602  * @brief Load energy boundaries from FITS file
603  *
604  * @param[in] filename FITS file name.
605  *
606  * Loads the energy boundaries from a FITS file.
607  *
608  * If no extension name is provided, the energy boundaries are loaded from
609  * the `EBOUNDS` extension.
610  ***************************************************************************/
611 void GEbounds::load(const GFilename& filename)
612 {
613  // Open FITS file
614  GFits fits(filename);
615 
616  // Get energy boundary table
617  const GFitsTable& table = *fits.table(filename.extname(gammalib::extname_ebounds));
618 
619  // Read energy boundaries from table
620  read(table);
621 
622  // Close FITS file
623  fits.close();
624 
625  // Return
626  return;
627 }
628 
629 
630 /***********************************************************************//**
631  * @brief Save energy boundaries into FITS file
632  *
633  * @param[in] filename FITS file name.
634  * @param[in] clobber Overwrite an existing energy boundaries extension?
635  * @param[in] unit Energy unit
636  *
637  * Saves energy boundaries into a FITS file. If a file with the given
638  * @p filename does not yet exist it will be created, otherwise the method
639  * opens the existing file. Energy boundaries can only be appended to an
640  * existing file if the @p clobber flag is set to "true" (otherwise an
641  * exception is thrown).
642  *
643  * The method will append a binary FITS table containing the energy
644  * boundaries to the FITS file. The extension name can be specified as part
645  * of the @p filename. For example the @p filename
646  *
647  * myfile.fits[ENERGY BOUNDARIES]
648  *
649  * will save the energy boundaries in the `ENERGY BOUNDARIES` extension of
650  * the "myfile.fits" file. If the extension exists already in the file it
651  * will be replaced, otherwise a new extension will be created. If no
652  * extension name is provided, the method will use `EBOUNDS` as the default
653  * extension name for energy boundaries.
654  ***************************************************************************/
655 void GEbounds::save(const GFilename& filename,
656  const bool& clobber,
657  const std::string& unit) const
658 {
659  // Open or create FITS file (without extension name since the requested
660  // extension may not yet exist in the file)
661  GFits fits(filename.url(), true);
662 
663  // Write energy boundaries to FITS file
664  write(fits, filename.extname(gammalib::extname_ebounds), unit);
665 
666  // Save to file
667  fits.save(clobber);
668 
669  // Return
670  return;
671 }
672 
673 
674 /***********************************************************************//**
675  * @brief Read energy boundaries from FITS table
676  *
677  * @param[in] table FITS table.
678  *
679  * Reads the energy boundaries from a FITS table. The method interprets the
680  * energy units provide in the FITS header. If no energy units are found it
681  * is assumed that the energies are stored in units of keV.
682  ***************************************************************************/
683 void GEbounds::read(const GFitsTable& table)
684 {
685  // Free members
686  free_members();
687 
688  // Initialise attributes
689  init_members();
690 
691  // Extract energy boundary information from FITS table
692  m_num = table.integer("NAXIS2");
693  if (m_num > 0) {
694 
695  // Allocate memory
696  m_min = new GEnergy[m_num];
697  m_max = new GEnergy[m_num];
698 
699  // Get units
700  std::string emin_unit = table["E_MIN"]->unit();
701  std::string emax_unit = table["E_MAX"]->unit();
702  if (emin_unit.empty()) {
703  emin_unit = "keV";
704  }
705  if (emax_unit.empty()) {
706  emax_unit = "keV";
707  }
708 
709  // Copy information
710  for (int i = 0; i < m_num; ++i) {
711  m_min[i](table["E_MIN"]->real(i), emin_unit);
712  m_max[i](table["E_MAX"]->real(i), emax_unit);
713  }
714 
715  // Set attributes
716  set_attributes();
717 
718  } // endif: there were channels to read
719 
720  // Return
721  return;
722 }
723 
724 
725 /***********************************************************************//**
726  * @brief Write energy boundaries into FITS object
727  *
728  * @param[in] fits FITS file.
729  * @param[in] extname Energy boundary extension name.
730  * @param[in] unit Energy units.
731  *
732  * Writes the energy boundaries into a FITS object. The @p unit parameter
733  * specifies in which unit the energies are written. By default, the energy
734  * units are keV.
735  *
736  * @todo Write header keywords.
737  ***************************************************************************/
739  const std::string& extname,
740  const std::string& unit) const
741 {
742  // Create energy boundary columns
743  GFitsTableDoubleCol cemin("E_MIN", m_num);
744  GFitsTableDoubleCol cemax("E_MAX", m_num);
745 
746  // Fill energy boundary columns
747  for (int i = 0; i < m_num; ++i) {
748  cemin(i) = m_min[i](unit);
749  cemax(i) = m_max[i](unit);
750  }
751 
752  // Set energy units
753  cemin.unit(unit);
754  cemax.unit(unit);
755 
756  // Create binary table
757  GFitsBinTable table(m_num);
758  table.append(cemin);
759  table.append(cemax);
760  table.extname(extname);
761 
762  // If the FITS object contains already an extension with the same
763  // name then remove now this extension
764  if (fits.contains(extname)) {
765  fits.remove(extname);
766  }
767 
768  // Append energy boundary table to FITS file
769  fits.append(table);
770 
771  // Return
772  return;
773 }
774 
775 
776 /***********************************************************************//**
777  * @brief Read energy boundaries from XML element
778  *
779  * @param[in] xml XML element.
780  *
781  * @exception GException::invalid_value
782  * Invalid XML format encountered.
783  *
784  * Read energy boundaries from an XML element. The format of the energy
785  * boundaries is
786  *
787  * <parameter name="EnergyBoundaries" emin="0.1" emax="10.0"/>
788  *
789  * The units of the @a emin and @a emax parameters are MeV.
790  ***************************************************************************/
791 void GEbounds::read(const GXmlElement& xml)
792 {
793  // Clear energy boundaries
794  clear();
795 
796  // Get energy boundaries parameter
797  const GXmlElement* par = gammalib::xml_get_par(G_READ_XML, xml,
798  "EnergyBoundaries");
799 
800  // Extract position attributes
801  if (par->has_attribute("emin") && par->has_attribute("emax")) {
802  double emin = gammalib::todouble(par->attribute("emin"));
803  double emax = gammalib::todouble(par->attribute("emax"));
804  append(GEnergy(emin, "MeV"), GEnergy(emax, "MeV"));
805  }
806  else {
807  std::string msg = "Attributes \"emin\" and/or \"emax\" not found"
808  " in XML parameter \"EnergyBoundaries\"."
809  " Please verify the XML format.";
811  }
812 
813  // Set attribues
814  set_attributes();
815 
816  // Return
817  return;
818 }
819 
820 
821 /***********************************************************************//**
822  * @brief Write energy boundaries into XML element
823  *
824  * @param[in] xml XML element.
825  *
826  * Writes energy boundaries into an XML element. The format of the energy
827  * boundaries is
828  *
829  * <parameter name="EnergyBoundaries" emin="0.1" emax="10.0"/>
830  *
831  * The units of the @a emin and @a emax parameters are MeV.
832  *
833  * This method does nothing if the energy boundaries are empty.
834  ***************************************************************************/
835 void GEbounds::write(GXmlElement& xml) const
836 {
837  // Continue only if there are energy boundaries
838  if (!is_empty()) {
839 
840  // Get parameter
842  "EnergyBoundaries");
843 
844  // Write attributes
845  par->attribute("emin", gammalib::str(emin().MeV()));
846  par->attribute("emax", gammalib::str(emax().MeV()));
847 
848  } // endif: energy boundaries were not empty
849 
850  // Return
851  return;
852 }
853 
854 
855 /***********************************************************************//**
856  * @brief Returns energy bin index for a given energy
857  *
858  * @param[in] eng Energy.
859  * @return Bin index.
860  *
861  * Returns the energy boundary bin index for a given energy. By convention,
862  * the limits for an energy bin are defined as
863  *
864  * min <= energy < max
865  *
866  * i.e. and energy equals to max falls above the largest energy.
867  *
868  * If the energy falls outside all boundaries, -1 is returned.
869  ***************************************************************************/
870 int GEbounds::index(const GEnergy& eng) const
871 {
872  // Initialise index with 'not found'
873  int index = -1;
874 
875  // Search all energy boundaries for containment
876  for (int i = 0; i < m_num; ++i) {
877  if (eng >= m_min[i] && eng < m_max[i]) {
878  index = i;
879  break;
880  }
881  }
882 
883  // Return index
884  return index;
885 }
886 
887 
888 /***********************************************************************//**
889  * @brief Set minimum energy for a given energy interval
890  *
891  * @param[in] index Energy interval index (0,...,size()-1).
892  * @param[in] energy Minimum energy of interval.
893  *
894  * @exception GException::out_of_range
895  * Specified index is out of range.
896  *
897  * Sets the minimum energy for the energy interval @p index.
898  ***************************************************************************/
899 void GEbounds::emin(const int& index, const GEnergy& energy)
900 {
901  #if defined(G_RANGE_CHECK)
902  // Throw an exception if index is outside valid range
903  if (index < 0 || index >= m_num) {
905  "Minimum energy of interval",
906  index, m_num);
907  }
908  #endif
909 
910  // Set minimum energy
911  m_min[index] = energy;
912 
913  // Set attributes
914  set_attributes();
915 
916  // Return
917  return;
918 }
919 
920 
921 /***********************************************************************//**
922  * @brief Set maximum energy for a given energy interval
923  *
924  * @param[in] index Energy interval index (0,...,size()-1).
925  * @param[in] energy Maximum energy of interval.
926  *
927  * @exception GException::out_of_range
928  * Specified index is out of range.
929  *
930  * Sets the maximum energy for the energy interval @p index.
931  ***************************************************************************/
932 void GEbounds::emax(const int& index, const GEnergy& energy)
933 {
934  #if defined(G_RANGE_CHECK)
935  // Throw an exception if index is outside valid range
936  if (index < 0 || index >= m_num) {
938  "Maximum energy of interval",
939  index, m_num);
940  }
941  #endif
942 
943  // Set maximum energy
944  m_max[index] = energy;
945 
946  // Set attributes
947  set_attributes();
948 
949  // Return
950  return;
951 }
952 
953 
954 /***********************************************************************//**
955  * @brief Returns minimum energy for a given energy interval
956  *
957  * @param[in] index Energy interval index (0,...,size()-1).
958  * @return Minimum energy of interval.
959  *
960  * @exception GException::out_of_range
961  * Specified index is out of range.
962  ***************************************************************************/
963 GEnergy GEbounds::emin(const int& index) const
964 {
965  #if defined(G_RANGE_CHECK)
966  // Throw an exception if index is outside valid range
967  if (index < 0 || index >= m_num) {
969  "Minimum energy of interval",
970  index, m_num);
971  }
972  #endif
973 
974  // Return minimum energy
975  return (m_min[index]);
976 }
977 
978 
979 /***********************************************************************//**
980  * @brief Returns maximum energy for a given energy interval
981  *
982  * @param[in] index Energy interval index (0,...,size()-1).
983  * @return Maximum energy of interval.
984  *
985  * @exception GException::out_of_range
986  * Specified index is out of range.
987  ***************************************************************************/
988 GEnergy GEbounds::emax(const int& index) const
989 {
990  #if defined(G_RANGE_CHECK)
991  // Throw an exception if index is outside valid range
992  if (index < 0 || index >= m_num) {
994  "Maximum energy of interval",
995  index, m_num);
996  }
997  #endif
998 
999  // Return maximum energy
1000  return (m_max[index]);
1001 }
1002 
1003 
1004 /***********************************************************************//**
1005  * @brief Returns mean energy for a given energy interval
1006  *
1007  * @param[in] index Energy interval index (0,...,size()-1).
1008  * @return Mean energy of interval.
1009  *
1010  * @exception GException::out_of_range
1011  * Specified index is out of range.
1012  *
1013  * Computes the mean energy
1014  * \f$0.5 * (E_{\rm min} + E_{\rm max})\f$
1015  * for the energy interval @p index.
1016  ***************************************************************************/
1017 GEnergy GEbounds::emean(const int& index) const
1018 {
1019  #if defined(G_RANGE_CHECK)
1020  // If index is outside boundary then throw an error
1021  if (index < 0 || index >= m_num) {
1022  throw GException::out_of_range(G_EMEAN, index, 0, m_num-1);
1023  }
1024  #endif
1025 
1026  // Compute mean energy
1027  GEnergy emean = 0.5 * (m_min[index] + m_max[index]);
1028 
1029  // Return
1030  return emean;
1031 }
1032 
1033 
1034 /***********************************************************************//**
1035  * @brief Returns logarithmic mean energy for a given energy interval
1036  *
1037  * @param[in] index Energy interval index (0,...,size()-1).
1038  * @return Logarithmic mean energy of interval.
1039  *
1040  * @exception GException::out_of_range
1041  * Specified index is out of range.
1042  *
1043  * Computes the logarithmic mean energy
1044  * \f$10^{0.5 * (\log E_{\rm min} + \log E_{\rm max})}\f$
1045  * for the energy interval @p index.
1046  ***************************************************************************/
1047 GEnergy GEbounds::elogmean(const int& index) const
1048 {
1049  #if defined(G_RANGE_CHECK)
1050  // If index is outside boundary then throw an error
1051  if (index < 0 || index >= m_num) {
1052  throw GException::out_of_range(G_ELOGMEAN, index, 0, m_num-1);
1053  }
1054  #endif
1055 
1056  // Compute logarithmic mean energy
1057  GEnergy elogmean;
1058  elogmean.MeV(std::sqrt(m_min[index].MeV() * m_max[index].MeV()));
1059 
1060  // Return
1061  return elogmean;
1062 }
1063 
1064 
1065 /***********************************************************************//**
1066  * @brief Returns energy interval width
1067  *
1068  * @param[in] index Energy interval index (0,...,size()-1).
1069  * @return Energy interval width.
1070  *
1071  * @exception GException::out_of_range
1072  * Specified index is out of range.
1073  ***************************************************************************/
1074 GEnergy GEbounds::ewidth(const int& index) const
1075 {
1076  #if defined(G_RANGE_CHECK)
1077  // If index is outside boundary then throw an error
1078  if (index < 0 || index >= m_num) {
1079  throw GException::out_of_range(G_EWIDTH, index, 0, m_num-1);
1080  }
1081  #endif
1082 
1083  // Return
1084  return (m_max[index]-m_min[index]);
1085 }
1086 
1087 
1088 
1089 
1090 /***********************************************************************//**
1091  * @brief Checks whether energy boundaries contain energy
1092  *
1093  * @param[in] eng Energy to be checked.
1094  * @return True if energy falls in at least one interval, false otherwise.
1095  *
1096  * Checks if the energy @p eng falls in at least one of the energy intervals.
1097  * The method exits when the first matching interval has been found.
1098  ***************************************************************************/
1099 bool GEbounds::contains(const GEnergy& eng) const
1100 {
1101  // Initialise test
1102  bool found = false;
1103 
1104  // Test all energy boundaries
1105  for (int i = 0; i < m_num; ++i) {
1106  if (eng >= m_min[i] && eng <= m_max[i]) {
1107  found = true;
1108  break;
1109  }
1110  }
1111 
1112  // Return result
1113  return found;
1114 }
1115 
1116 
1117 /***********************************************************************//**
1118  * @brief Checks whether energy boundaries contain and energy bin
1119  *
1120  * @param[in] emin Minimum energy of bin to be checked.
1121  * @param[in] emax Maximum energy of bin to be checked.
1122  * @return True if energy bin [emin, emax] is fully contained inside energy
1123  * the boundaries, false otherwiese
1124  *
1125  * Checks if the energy interval [@p emin, @p emax ] falls is fully contained
1126  * by energy boundaries
1127  *
1128  * @todo This method is so far only correct for contiguous energy boundaries.
1129  ***************************************************************************/
1130 bool GEbounds::contains(const GEnergy& emin, const GEnergy& emax) const
1131 {
1132  // Initialise result
1133  bool contained = false;
1134 
1135  // Check if emin, emax is fully contained within this ebounds
1136  if (emin >= m_emin && emax <= m_emax) {
1137  contained = true;
1138  }
1139 
1140  // Return result
1141  return contained;
1142 }
1143 
1144 
1145 /***********************************************************************//**
1146  * @brief Print energy boundaries
1147  *
1148  * @param[in] chatter Chattiness (defaults to NORMAL).
1149  * @return String containing energy boundary information.
1150  ***************************************************************************/
1151 std::string GEbounds::print(const GChatter& chatter) const
1152 {
1153  // Initialise result string
1154  std::string result;
1155 
1156  // Continue only if chatter is not silent
1157  if (chatter != SILENT) {
1158 
1159  // Append Header
1160  result.append("=== GEbounds ===");
1161 
1162  // Append information
1163  result.append("\n"+gammalib::parformat("Number of intervals"));
1164  result.append(gammalib::str(size()));
1165  result.append("\n"+gammalib::parformat("Energy range"));
1166  result.append(emin().print());
1167  result.append(" - ");
1168  result.append(emax().print());
1169 
1170  // If there are multiple energy bins then append them
1171  if (chatter >= EXPLICIT) {
1172  if (size() > 1) {
1173  for (int i = 0; i < size(); ++i) {
1174  result.append("\n");
1175  result.append(gammalib::parformat("Energy interval "+
1176  gammalib::str(i)));
1177  result.append(emin(i).print());
1178  result.append(" - ");
1179  result.append(emax(i).print());
1180  }
1181  }
1182  } // endif: chatter was less than explicit
1183 
1184  } // endif: chatter was not silent
1185 
1186  // Return
1187  return result;
1188 }
1189 
1190 
1191 /*==========================================================================
1192  = =
1193  = Private methods =
1194  = =
1195  ==========================================================================*/
1196 
1197 /***********************************************************************//**
1198  * @brief Initialise class members
1199  ***************************************************************************/
1201 {
1202  // Initialise members
1203  m_num = 0;
1204  m_emin.clear();
1205  m_emax.clear();
1206  m_min = NULL;
1207  m_max = NULL;
1208 
1209  // Return
1210  return;
1211 }
1212 
1213 
1214 /***********************************************************************//**
1215  * @brief Copy class members
1216  *
1217  * @param[in] ebds Energy boundaries.
1218  ***************************************************************************/
1220 {
1221  // Copy attributes
1222  m_num = ebds.m_num;
1223  m_emin = ebds.m_emin;
1224  m_emax = ebds.m_emax;
1225 
1226  // Copy arrays
1227  if (m_num > 0) {
1228  m_min = new GEnergy[m_num];
1229  m_max = new GEnergy[m_num];
1230  for (int i = 0; i < m_num; ++i) {
1231  m_min[i] = ebds.m_min[i];
1232  m_max[i] = ebds.m_max[i];
1233  }
1234  }
1235 
1236  // Return
1237  return;
1238 }
1239 
1240 
1241 /***********************************************************************//**
1242  * @brief Delete class members
1243  ***************************************************************************/
1245 {
1246  // Free memory
1247  if (m_min != NULL) delete [] m_min;
1248  if (m_max != NULL) delete [] m_max;
1249 
1250  // Signal free pointers
1251  m_min = NULL;
1252  m_max = NULL;
1253 
1254  // Return
1255  return;
1256 }
1257 
1258 
1259 /***********************************************************************//**
1260  * @brief Set class attributes
1261  *
1262  * Determines the minimum and maximum energy from all intervals. If no
1263  * interval is present the minimum and maximum energies are cleared.
1264  ***************************************************************************/
1266 {
1267  // If there are intervals then determine the minimum and maximum
1268  // energy from these intervals ...
1269  if (m_num > 0) {
1270  m_emin = m_min[0];
1271  m_emax = m_max[0];
1272  for (int i = 1; i < m_num; ++i) {
1273  if (m_min[i] < m_emin) m_emin = m_min[i];
1274  if (m_max[i] > m_emax) m_emax = m_max[i];
1275  }
1276  }
1277 
1278  // ... otherwise clear the minimum and maximum energy
1279  else {
1280  m_emin.clear();
1281  m_emax.clear();
1282  }
1283 
1284  // Return
1285  return;
1286 }
1287 
1288 
1289 /***********************************************************************//**
1290  * @brief Insert energy interval
1291  *
1292  * @param[in] index Index after with interval is inserted.
1293  * @param[in] emin Minimum energy of interval.
1294  * @param[in] emax Maximum energy of interval.
1295  *
1296  * @exception GException::invalid_argument
1297  * Minimum energy larger than maximum energy
1298  *
1299  * Inserts an energy interval after the specified @p index in the energy
1300  * boundaries. The method does not reorder the intervals by energy, instead
1301  * the client needs to determine the approriate @p index.
1302  *
1303  * Invalid parameters do not produce any exception, but are handled
1304  * transparently. If the interval is invalid (i.e. @p emin > @p emax) an
1305  * exception is thrown. If the @p index is out of the valid range, the
1306  * index will be adjusted to either the first or the last element.
1307  ***************************************************************************/
1308 void GEbounds::insert_eng(const int& index,
1309  const GEnergy& emin,
1310  const GEnergy& emax)
1311 {
1312  // Throw an exception if energy interval is invalid
1313  if (emin > emax) {
1314  std::string msg = "Invalid energy interval specified. Minimum"
1315  " energy "+emin.print(NORMAL)+" can not be"
1316  " larger than maximum energy "+
1317  emax.print(NORMAL)+".";
1319  }
1320 
1321  // Set index
1322  int inx = index;
1323 
1324  // If inx is out of range then adjust it
1325  if (inx < 0) inx = 0;
1326  if (inx > m_num) inx = m_num;
1327 
1328  // Allocate new intervals
1329  int num = m_num+1;
1330  GEnergy* min = new GEnergy[num];
1331  GEnergy* max = new GEnergy[num];
1332 
1333  // Copy intervals before index to be inserted
1334  for (int i = 0; i < inx; ++i) {
1335  min[i] = m_min[i];
1336  max[i] = m_max[i];
1337  }
1338 
1339  // Insert interval
1340  min[inx] = emin;
1341  max[inx] = emax;
1342 
1343  // Copy intervals after index to be inserted
1344  for (int i = inx+1; i < num; ++i) {
1345  min[i] = m_min[i-1];
1346  max[i] = m_max[i-1];
1347  }
1348 
1349  // Free memory
1350  if (m_min != NULL) delete [] m_min;
1351  if (m_max != NULL) delete [] m_max;
1352 
1353  // Set new memory
1354  m_min = min;
1355  m_max = max;
1356 
1357  // Set number of elements
1358  m_num = num;
1359 
1360  // Set attributes
1361  set_attributes();
1362 
1363  // Return
1364  return;
1365 }
1366 
1367 
1368 /*==========================================================================
1369  = =
1370  = Friends =
1371  = =
1372  ==========================================================================*/
1373 
1374 /***********************************************************************//**
1375  * @brief Energy boundaries equality operator friend
1376  *
1377  * @param[in] a First energy boundaries.
1378  * @param[in] b Second energy boundaries.
1379  * @return True if both energy boundaries are identical.
1380  ***************************************************************************/
1381 bool operator==(const GEbounds& a, const GEbounds& b)
1382 {
1383  // Initialise identify flag
1384  bool identity = true;
1385 
1386  // Check that both energy boundaries have the same size
1387  if (a.size() != b.size()) {
1388  identity = false;
1389  }
1390 
1391  // Check all energy boundaries
1392  else {
1393  for (int i = 0; i < a.size(); ++i) {
1394  if ((a.emin(i) != b.emin(i)) || (a.emax(i) != b.emax(i))) {
1395  identity = false;
1396  break;
1397  }
1398  }
1399  }
1400 
1401  // Return identity flag
1402  return identity;
1403 }
GEnergy * m_max
Array of interval maximum energies.
Definition: GEbounds.hpp:140
void unit(const std::string &unit)
Set column unit.
void init_members(void)
Initialise class members.
Definition: GEbounds.cpp:1200
FITS table double column class interface definition.
#define G_WRITE_XML
Definition: GEbounds.cpp:45
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
bool contains(const int &extno) const
Check if HDU exists in FITS file.
Definition: GFits.hpp:282
void copy_members(const GEbounds &ebds)
Copy class members.
Definition: GEbounds.cpp:1219
GEnergy m_emin
Minimum energy of all intervals.
Definition: GEbounds.hpp:137
XML element node class interface definition.
GEbounds(void)
Constructor.
Definition: GEbounds.cpp:72
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:162
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition: GEbounds.cpp:300
GEbounds * clone(void) const
Clone energy boundaries.
Definition: GEbounds.cpp:285
#define G_READ_XML
Definition: GEbounds.cpp:44
bool contains(const GEnergy &eng) const
Checks whether energy boundaries contain energy.
Definition: GEbounds.cpp:1099
void set_attributes(void)
Set class attributes.
Definition: GEbounds.cpp:1265
XML element node class.
Definition: GXmlElement.hpp:47
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition: GFitsTable.hpp:147
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
void read(const GFitsTable &table)
Read energy boundaries from FITS table.
Definition: GEbounds.cpp:683
#define G_EMAX_SET
Definition: GEbounds.cpp:48
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
GEnergy emean(const int &index) const
Returns mean energy for a given energy interval.
Definition: GEbounds.cpp:1017
FITS file class interface definition.
double min(const GVector &vector)
Computes minimum vector element.
Definition: GVector.cpp:850
GEnergy ewidth(const int &index) const
Returns energy interval width.
Definition: GEbounds.cpp:1074
void write(GFits &file, const std::string &extname=gammalib::extname_ebounds, const std::string &unit="keV") const
Write energy boundaries into FITS object.
Definition: GEbounds.cpp:738
Energy container class.
Definition: GEnergies.hpp:60
int index(const GEnergy &eng) const
Returns energy bin index for a given energy.
Definition: GEbounds.cpp:870
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
int m_num
Number of energy boundaries.
Definition: GEbounds.hpp:136
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1275
void extend(const GEbounds &ebds)
Append energy boundaries.
Definition: GEbounds.cpp:474
#define G_EMEAN
Definition: GEbounds.cpp:51
#define G_EWIDTH
Definition: GEbounds.cpp:53
void insert_eng(const int &index, const GEnergy &emin, const GEnergy &emax)
Insert energy interval.
Definition: GEbounds.cpp:1308
void insert(const GEnergy &emin, const GEnergy &emax)
Insert energy interval.
Definition: GEbounds.cpp:321
Energy boundaries container class.
Definition: GEbounds.hpp:60
#define G_INSERT_ENG
Definition: GEbounds.cpp:54
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:848
const GXmlAttribute * attribute(const int &index) const
Return attribute.
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
void set(const GEnergies &energies)
Set energy boundaries from energy container.
Definition: GEbounds.cpp:540
Filename class.
Definition: GFilename.hpp:62
Energy container class definition.
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:186
bool has_attribute(const std::string &name) const
Check if element has a given attribute.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:170
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1527
#define G_REMOVE
Definition: GEbounds.cpp:46
bool is_empty(void) const
Signal if there are no energy boundaries.
Definition: GEbounds.hpp:174
const std::string extname_ebounds
Definition: GEbounds.hpp:44
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
std::string print(const GChatter &chatter=NORMAL) const
Print energy boundaries.
Definition: GEbounds.cpp:1151
GChatter
Definition: GTypemaps.hpp:33
GEnergy elogmean(const int &index) const
Returns logarithmic mean energy for a given energy interval.
Definition: GEbounds.cpp:1047
int integer(const std::string &keyname) const
Return card value as integer.
Definition: GFitsHDU.hpp:436
void reserve(const int &num)
Reserve space for energy intervals.
Definition: GEbounds.cpp:460
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:267
double max(const GVector &vector)
Computes maximum vector element.
Definition: GVector.cpp:879
#define G_EMIN_GET
Definition: GEbounds.cpp:49
GEnergy * m_min
Array of interval minimum energies.
Definition: GEbounds.hpp:139
void free_members(void)
Delete class members.
Definition: GEbounds.cpp:1244
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
GEnergy m_emax
Maximum energy of all intervals.
Definition: GEbounds.hpp:138
#define G_EMIN_SET
Definition: GEbounds.cpp:47
virtual ~GEbounds(void)
Destructor.
Definition: GEbounds.cpp:215
FITS binary table class.
Energy boundaries class interface definition.
void remove(const int &index)
Remove energy interval.
Definition: GEbounds.cpp:426
Exception handler interface definition.
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition: GFits.cpp:665
FITS binary table class definition.
#define G_ELOGMEAN
Definition: GEbounds.cpp:52
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:198
bool operator==(const GEnergy &a, const GEnergy &b)
Energy equality operator friend.
Definition: GEnergy.hpp:297
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1033
#define G_EMAX_GET
Definition: GEbounds.cpp:50
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
void save(const GFilename &filename, const bool &clobber=false, const std::string &unit="keV") const
Save energy boundaries into FITS file.
Definition: GEbounds.cpp:655
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1579
GEbounds & operator=(const GEbounds &ebds)
Assignment operator.
Definition: GEbounds.cpp:237
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
FITS table double column.
void merge(void)
Merge all overlapping or connecting successive energy intervals.
Definition: GEbounds.cpp:348
Filename class interface definition.
double todouble(const std::string &arg)
Convert string into double precision value.
Definition: GTools.cpp:816
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:415
FITS table abstract base class interface definition.
void load(const GFilename &filename)
Load energy boundaries from FITS file.
Definition: GEbounds.cpp:611