GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GPha.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GPha.cpp - XSPEC Pulse Height Analyzer class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2013-2018 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 GPha.cpp
23  * @brief XSPEC Pulse Height Analyzer class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GPha.hpp"
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GEnergy.hpp"
35 #include "GFits.hpp"
36 #include "GFitsTable.hpp"
37 #include "GFitsBinTable.hpp"
38 #include "GFitsTableShortCol.hpp"
39 #include "GFitsTableFloatCol.hpp"
40 #include "GNdarray.hpp"
41 
42 /* __ Method name definitions ____________________________________________ */
43 #define G_OPERATOR_PLUS "GPha::operator+=(GPha&)"
44 #define G_OPERATOR_MINUS "GPha::operator-=(GPha&)"
45 #define G_OPERATOR "GPha::operator[](std::string&)"
46 #define G_AT1 "GPha::at(int&)"
47 #define G_AT2 "GPha::at(int&, int&)"
48 #define G_APPEND "GPha::append(std::string&, std::vector<double>&)"
49 #define G_AREASCAL_SET "GPha::areascal(int&, double&)"
50 #define G_AREASCAL_GET "GPha::areascal(int&)"
51 #define G_BACKSCAL_SET "GPha::backscal(int&, double&)"
52 #define G_BACKSCAL_GET "GPha::backscal(int&)"
53 #define G_READ "GPha::read(GFitsTable*)"
54 
55 /* __ Macros _____________________________________________________________ */
56 
57 /* __ Coding definitions _________________________________________________ */
58 
59 /* __ Debug definitions __________________________________________________ */
60 
61 
62 /*==========================================================================
63  = =
64  = Constructors/destructors =
65  = =
66  ==========================================================================*/
67 
68 /***********************************************************************//**
69  * @brief Void constructor
70  ***************************************************************************/
72 {
73  // Initialise members
74  init_members();
75 
76  // Return
77  return;
78 }
79 
80 
81 /***********************************************************************//**
82  * @brief File constructor
83  *
84  * @param[in] filename File name.
85  ***************************************************************************/
86 GPha::GPha(const GFilename& filename)
87 {
88  // Initialise members
89  init_members();
90 
91  // Load PHA file
92  load(filename);
93 
94  // Return
95  return;
96 }
97 
98 
99 /***********************************************************************//**
100  * @brief Energy boundary constructor
101  *
102  * @param[in] ebds Energy boundaries.
103  ***************************************************************************/
104 GPha::GPha(const GEbounds& ebds)
105 {
106  // Initialise members
107  init_members();
108 
109  // Set energy boundaries
110  m_ebounds = ebds;
111 
112  // Initialize spectrum
113  alloc(ebds.size());
114 
115  // Return
116  return;
117 }
118 
119 
120 /***********************************************************************//**
121  * @brief Energy bins constructor
122  *
123  * @param[in] bins Number of energy bins.
124  ***************************************************************************/
125 GPha::GPha(const int& bins)
126 {
127  // Initialise members
128  init_members();
129 
130  // Initialize spectrum
131  alloc(bins);
132 
133  // Return
134  return;
135 }
136 
137 
138 /***********************************************************************//**
139  * @brief Copy constructor
140  *
141  * @param[in] pha Pulse Height Analyzer spectrum.
142  ***************************************************************************/
143 GPha::GPha(const GPha& pha)
144 {
145  // Initialise members
146  init_members();
147 
148  // Copy members
149  copy_members(pha);
150 
151  // Return
152  return;
153 }
154 
155 
156 /***********************************************************************//**
157  * @brief Destructor
158  ***************************************************************************/
160 {
161  // Free members
162  free_members();
163 
164  // Return
165  return;
166 }
167 
168 
169 /*==========================================================================
170  = =
171  = Operators =
172  = =
173  ==========================================================================*/
174 
175 /***********************************************************************//**
176  * @brief Assignment operator
177  *
178  * @param[in] pha Pulse Height Analyzer spectrum.
179  * @return Pulse Height Analyzer spectrum.
180  ***************************************************************************/
182 {
183  // Execute only if object is not identical
184  if (this != &pha) {
185 
186  // Free members
187  free_members();
188 
189  // Initialise members
190  init_members();
191 
192  // Copy members
193  copy_members(pha);
194 
195  } // endif: object was not identical
196 
197  // Return
198  return *this;
199 }
200 
201 
202 /***********************************************************************//**
203  * @brief Add spectrum
204  *
205  * @param[in] pha Pulse Height Analyzer spectrum.
206  * @return Sum of Pulse Height Analyzer spectra.
207  *
208  * @exception GException::invalid_value
209  * Incompatible spectrum.
210  *
211  * Adds the counts of a spectrum to the counts of the current spectrum. The
212  * operator also adds the exposure times of the spectrum to the current
213  * exposure time. The area scaling factor \f$\alpha\f$ is recomputed for
214  * each spectral bin using
215  *
216  * \f[
217  * \alpha = \frac{N_1 + N_2}{\frac{N_1}{\alpha_1} + \frac{N_2}{\alpha_2}}
218  * \f]
219  *
220  * where
221  * \f$N_1\f$ and \f$N_2\f$ are the number of events in the bin for spectrum
222  * 1 and 2, respectively, and
223  * \f$\alpha_1\f$ and \f$\alpha_2\f$ are the corresponding area scaling
224  * factors.
225  *
226  * The background scaling factor is not altered.
227  *
228  * The operator only works if the provide specturm has the same energy
229  * binning than the current spectrum.
230  ***************************************************************************/
232 {
233  // Throw an exception if the spectra are not compatible
234  if (this->ebounds() != pha.ebounds()) {
235  std::string msg = "Incompatible energy binning of Pulse Height "
236  "Analyzer spectrum.";
238  }
239 
240  // Add spectra
241  for (int i = 0; i < this->size(); ++i) {
242 
243  // Update the area scaling
244  double n1 = m_counts[i];
245  double n2 = pha.m_counts[i];
246  double a1 = m_areascal[i];
247  double a2 = pha.m_areascal[i];
248  double f1 = (a1 > 0.0) ? n1/a1 : 0.0;
249  double f2 = (a2 > 0.0) ? n2/a2 : 0.0;
250  double f = f1 + f2;
251  m_areascal[i] = (f > 0.0) ? (n1+n2)/f : 1.0;
252 
253  // Add counts
254  m_counts[i] += pha.m_counts[i];
255  }
256 
257  // Add attributes
258  m_underflow += pha.m_underflow;
259  m_overflow += pha.m_overflow;
260  m_outflow += pha.m_outflow;
261  m_exposure += pha.m_exposure;
262 
263  // Return
264  return *this;
265 }
266 
267 
268 /***********************************************************************//**
269  * @brief Subtract spectrum
270  *
271  * @param[in] pha Pulse Height Analyzer spectrum.
272  * @return Difference of Pulse Height Analyzer spectra.
273  *
274  * @exception GException::invalid_value
275  * Incompatible spectrum.
276  *
277  * Subtracts the counts of a spectrum from the counts of the current
278  * spectrum. The operator also subtracts the exposure times of the spectrum
279  * from the current exposure time. The area scaling factor \f$\alpha\f$ is
280  * recomputed for each spectral bin using
281  *
282  * \f[
283  * \alpha = \frac{N_1 - N_2}{\frac{N_1}{\alpha_1} - \frac{N_2}{\alpha_2}}
284  * \f]
285  *
286  * where
287  * \f$N_1\f$ and \f$N_2\f$ are the number of events in the bin for spectrum
288  * 1 and 2, respectively, and
289  * \f$\alpha_1\f$ and \f$\alpha_2\f$ are the corresponding area scaling
290  * factors.
291  *
292  * The background scaling factor is not altered.
293  *
294  * The operator only works if the provide specturm has the same energy
295  * binning than the current spectrum.
296  ***************************************************************************/
298 {
299  // Throw an exception if the spectra are not compatible
300  if (this->ebounds() != pha.ebounds()) {
301  std::string msg = "Incompatible energy binning of Pulse Height "
302  "Analyzer spectrum.";
304  }
305 
306  // Subtract spectra
307  for (int i = 0; i < this->size(); ++i) {
308 
309  // Update the area scaling
310  double n1 = m_counts[i];
311  double n2 = pha.m_counts[i];
312  double a1 = m_areascal[i];
313  double a2 = pha.m_areascal[i];
314  double f1 = (a1 > 0.0) ? n1/a1 : 0.0;
315  double f2 = (a2 > 0.0) ? n2/a2 : 0.0;
316  double f = f1 - f2;
317  m_areascal[i] = (f > 0.0) ? (n1-n2)/f : 1.0;
318 
319  // Subtract counts
320  m_counts[i] -= pha.m_counts[i];
321  }
322 
323  // Subtract attributes
324  m_underflow -= pha.m_underflow;
325  m_overflow -= pha.m_overflow;
326  m_outflow -= pha.m_outflow;
327  m_exposure -= pha.m_exposure;
328 
329  // Return
330  return *this;
331 }
332 
333 
334 /***********************************************************************//**
335  * @brief Scale spectrum values
336  *
337  * @param[in] scale Scale factor.
338  * @return Scaled Pulse Height Analyzer spectrum.
339  *
340  * Multiplies the counts of a spectrum with a scale factor. The exposure
341  * time and area and background scale factors are not altered.
342  ***************************************************************************/
343 GPha& GPha::operator*=(const double& scale)
344 {
345  // Scale spectrums
346  for (int i = 0; i < this->size(); ++i) {
347  m_counts[i] *= scale;
348  }
349 
350  // Scale attributes
351  m_underflow *= scale;
352  m_overflow *= scale;
353  m_outflow *= scale;
354 
355  // Return
356  return *this;
357 }
358 
359 
360 /***********************************************************************//**
361  * @brief Divide spectrum values
362  *
363  * @param[in] scale Division factor.
364  * @return Divided Pulse Height Analyzer spectrum.
365  *
366  * Divides the counts of a spectrum by a scale factor. The exposure time
367  * and the area and background scale factors are not altered.
368  ***************************************************************************/
369 GPha& GPha::operator/=(const double& scale)
370 {
371  // Scale spectrums
372  for (int i = 0; i < this->size(); ++i) {
373  m_counts[i] /= scale;
374  }
375 
376  // Scale attributes
377  m_underflow /= scale;
378  m_overflow /= scale;
379  m_outflow /= scale;
380 
381  // Return
382  return *this;
383 }
384 
385 
386 /***********************************************************************//**
387  * @brief Return additional vector column
388  *
389  * @param[in] colname Vector column name.
390  * @return Vector column.
391  *
392  * Returns reference to additional vector column.
393  ***************************************************************************/
394 std::vector<double>& GPha::operator[](const std::string& colname)
395 {
396  // Determine index of additional column (-1 if not found)
397  int index = column_index(colname);
398 
399  // Throw exception if index not found
400  if (index == -1) {
401  std::string msg = "Could not find additional column with name \""+
402  colname+"\".";
404  }
405 
406  // Return reference to vector column
407  return (m_coldata[index]);
408 }
409 
410 
411 /***********************************************************************//**
412  * @brief Return additional vector column (const version)
413  *
414  * @param[in] colname Vector column name.
415  * @return Vector column.
416  *
417  * Returns reference to additional vector column.
418  ***************************************************************************/
419 const std::vector<double>& GPha::operator[](const std::string& colname) const
420 {
421  // Determine index of additional column (-1 if not found)
422  int index = column_index(colname);
423 
424  // Throw exception if index not found
425  if (index == -1) {
426  std::string msg = "Could not find additional column with name \""+
427  colname+"\".";
429  }
430 
431  // Return reference to vector column
432  return (m_coldata[index]);
433 }
434 
435 
436 /*==========================================================================
437  = =
438  = Public methods =
439  = =
440  ==========================================================================*/
441 
442 /***********************************************************************//**
443  * @brief Clear object.
444  *
445  * Reset object to a clean initial state.
446  ***************************************************************************/
447 void GPha::clear(void)
448 {
449  // Free memory
450  free_members();
451 
452  // Initialise members
453  init_members();
454 
455  // Return
456  return;
457 }
458 
459 
460 /***********************************************************************//**
461  * @brief Clone object
462  *
463  * @return Pulse Height Analyzer spectrum.
464  ***************************************************************************/
465 GPha* GPha::clone(void) const
466 {
467  // Clone client
468  return new GPha(*this);
469 }
470 
471 
472 /***********************************************************************//**
473  * @brief Return content of spectral bin
474  *
475  * @param[in] index Bin index [0,...,size()-1].
476  *
477  * @exception GException::out_of_range
478  * Spectral bin index is out of range.
479  *
480  * Returns reference to content of spectral bin with specified @p index.
481  ***************************************************************************/
482 double& GPha::at(const int& index)
483 {
484  // Raise exception if index is out of range
485  if (index < 0 || index >= size()) {
486  throw GException::out_of_range(G_AT1, "Spectral bin index", index, size());
487  }
488 
489  // Return reference
490  return (m_counts[index]);
491 }
492 
493 
494 /***********************************************************************//**
495  * @brief Return content of spectral bin (const version)
496  *
497  * @param[in] index Bin index [0,...,size()-1].
498  *
499  * @exception GException::out_of_range
500  * Spectral bin index is out of range.
501  *
502  * Returns reference to content of spectral bin with specified @p index.
503  ***************************************************************************/
504 const double& GPha::at(const int& index) const
505 {
506  // Raise exception if index is out of range
507  if (index < 0 || index >= size()) {
508  throw GException::out_of_range(G_AT1, "Spectral bin index", index, size());
509  }
510 
511  // Return reference
512  return (m_counts[index]);
513 }
514 
515 
516 /***********************************************************************//**
517  * @brief Return content of additional columns
518  *
519  * @param[in] index Bin index [0,...,size()-1].
520  * @param[in] col Columns index [0,...,columns()-1].
521  *
522  * @exception GException::out_of_range
523  * Spectral bin or column index is out of range.
524  *
525  * Returns reference to content of additional columns.
526  ***************************************************************************/
527 double& GPha::at(const int& index, const int& col)
528 {
529  // Throw an exception if bin or column index is out of range
530  if (index < 0 || index >= size()) {
531  throw GException::out_of_range(G_AT2, "Bin index", index, size());
532  }
533  if (col < 0 || col >= columns()) {
534  throw GException::out_of_range(G_AT2, "Column index", col, columns());
535  }
536 
537  // Return reference
538  return (m_coldata[col][index]);
539 }
540 
541 
542 /***********************************************************************//**
543  * @brief Return content of additional columns (const version)
544  *
545  * @param[in] index Bin index [0,...,size()-1].
546  * @param[in] col Columns index [0,...,columns()-1].
547  *
548  * @exception GException::out_of_range
549  * Spectral bin or column index is out of range.
550  *
551  * Returns reference to content of additional columns.
552  ***************************************************************************/
553 const double& GPha::at(const int& index, const int& col) const
554 {
555  // Throw an exception if bin or column index is out of range
556  if (index < 0 || index >= size()) {
557  throw GException::out_of_range(G_AT2, "Bin index", index, size());
558  }
559  if (col < 0 || col >= columns()) {
560  throw GException::out_of_range(G_AT2, "Column index", col, columns());
561  }
562 
563  // Return reference
564  return (m_coldata[col][index]);
565 }
566 
567 
568 /***********************************************************************//**
569  * @brief Append additional column to spectrum
570  *
571  * @param[in] name Additional column name.
572  * @param[in] column Additional column data.
573  ***************************************************************************/
574 void GPha::append(const std::string& name, const std::vector<double>& column)
575 {
576  // Throw an exception if the number of elements in the column does not
577  // correspond to the size of the spectrum
578  if (column.size() != size()) {
579  std::string msg = "Size of column "+gammalib::str(column.size())+
580  " is incompatible with size of spectrum "+
581  gammalib::str(size())+".";
583  }
584 
585  // Append column name
586  m_colnames.push_back(name);
587 
588  // Append column data
589  m_coldata.push_back(column);
590 
591  // Return
592  return;
593 }
594 
595 
596 /***********************************************************************//**
597  * @brief Set area scaling factor
598  *
599  * @param[in] index Bin index [0,...,size()-1].
600  * @param[in] areascal Area scaling factor.
601  *
602  * @exception GException::out_of_range
603  * Spectral bin index is out of range.
604  *
605  * Set area scaling for spectral bin with specified @p index.
606  ***************************************************************************/
607 void GPha::areascal(const int& index, const double& areascal)
608 {
609  // Raise exception if index is out of range
610  if (index < 0 || index >= size()) {
611  throw GException::out_of_range(G_AREASCAL_SET, "Spectral bin index",
612  index, size());
613  }
614 
615  // Set area scaling factor
616  m_areascal[index] = areascal;
617 
618  // Return
619  return;
620 }
621 
622 
623 /***********************************************************************//**
624  * @brief Return area scaling factor
625  *
626  * @param[in] index Bin index [0,...,size()-1].
627  *
628  * @exception GException::out_of_range
629  * Spectral bin index is out of range.
630  *
631  * Returns reference to area scaling for spectral bin with specified
632  * @p index.
633  ***************************************************************************/
634 const double& GPha::areascal(const int& index) const
635 {
636  // Raise exception if index is out of range
637  if (index < 0 || index >= size()) {
638  throw GException::out_of_range(G_AREASCAL_GET, "Spectral bin index",
639  index, size());
640  }
641 
642  // Return reference
643  return (m_areascal[index]);
644 }
645 
646 
647 /***********************************************************************//**
648  * @brief Set background scaling factor
649  *
650  * @param[in] index Bin index [0,...,size()-1].
651  * @param[in] backscal Background scaling factor.
652  *
653  * @exception GException::out_of_range
654  * Spectral bin index is out of range.
655  *
656  * Set background scaling for spectral bin with specified @p index.
657  ***************************************************************************/
658 void GPha::backscal(const int& index, const double& backscal)
659 {
660  // Raise exception if index is out of range
661  if (index < 0 || index >= size()) {
662  throw GException::out_of_range(G_BACKSCAL_SET, "Spectral bin index",
663  index, size());
664  }
665 
666  // Set background scaling factor
667  m_backscal[index] = backscal;
668 
669  // Return
670  return;
671 }
672 
673 
674 /***********************************************************************//**
675  * @brief Return background scaling factor
676  *
677  * @param[in] index Bin index [0,...,size()-1].
678  *
679  * @exception GException::out_of_range
680  * Spectral bin index is out of range.
681  *
682  * Returns reference to background scaling for spectral bin with specified
683  * @p index.
684  ***************************************************************************/
685 const double& GPha::backscal(const int& index) const
686 {
687  // Raise exception if index is out of range
688  if (index < 0 || index >= size()) {
689  throw GException::out_of_range(G_BACKSCAL_GET, "Spectral bin index",
690  index, size());
691  }
692 
693  // Return reference
694  return (m_backscal[index]);
695 }
696 
697 
698 /***********************************************************************//**
699  * @brief Number of counts in spectrum
700  *
701  * Returns the number of counts in the spectrum.
702  ***************************************************************************/
703 double GPha::counts(void) const
704 {
705  // Initialise counts
706  double counts = 0.0;
707 
708  // Compute content
709  for (int i = 0; i < m_counts.size(); ++i) {
710  counts += m_counts[i];
711  }
712 
713  // Return counts
714  return counts;
715 }
716 
717 /***********************************************************************//**
718  * @brief Get number of counts in spectrum as GNdarray
719  *
720  * @return GNdarray of number of counts in spectrum.
721  *
722  * Returns the number of counts in the spectrum as GNdarray.
723  ***************************************************************************/
725 {
726  // Initialise array
727  int size = m_counts.size();
728  GNdarray counts = GNdarray(size);
729 
730  // Compute content
731  for (int i = 0; i < size; ++i) {
732  counts(i) = m_counts[i];
733  }
734 
735  // Return counts
736  return counts;
737 }
738 
739 
740 /***********************************************************************//**
741  * @brief Get background scaling factors as GNdarray
742  *
743  * @return GNdarray of background scaling factors.
744  *
745  * Returns the background scaling factors as GNdarray.
746  ***************************************************************************/
748 {
749  // Initialise array
750  int size = m_backscal.size();
751  GNdarray alpha = GNdarray(size);
752 
753  // Compute content
754  for (int i = 0; i < size; ++i) {
755  alpha(i) = m_backscal[i];
756  }
757 
758  // Return background scaling factors
759  return alpha;
760 }
761 
762 
763 /***********************************************************************//**
764  * @brief Fill spectrum with a value.
765  *
766  * @param[in] energy Energy.
767  * @param[in] value Fill value (defaults to 1.0).
768  *
769  * Fills the specified @p value at a given @p energy in the spectrum.
770  ***************************************************************************/
771 void GPha::fill(const GEnergy& energy, const double& value)
772 {
773  // Get index
774  int index = m_ebounds.index(energy);
775 
776  // If index is not valid, then check for underflow, overflow or outflow
777  if (index == -1) {
778  if (energy < m_ebounds.emin()) {
779  m_underflow += value;
780  }
781  else if (energy >= m_ebounds.emax()) {
782  m_overflow += value;
783  }
784  else {
785  m_outflow += value;
786  }
787  }
788 
789  // ... otherwise fill the histogram
790  else {
791  m_counts[index] += value;
792  }
793 
794  // Return
795  return;
796 }
797 
798 
799 /***********************************************************************//**
800  * @brief Load Pulse Height Analyzer spectrum
801  *
802  * @param[in] filename File name.
803  *
804  * Loads the Pulse Height Analyzer spectrum from a FITS file.
805  ***************************************************************************/
806 void GPha::load(const GFilename& filename)
807 {
808  // Clear spectrum
809  clear();
810 
811  // Open FITS file (without extension name as the user is not allowed
812  // to modify the extension names)
813  GFits fits(filename.url());
814 
815  // Read PHA data
816  read(fits);
817 
818  // Close FITS file
819  fits.close();
820 
821  // Store filename
822  m_filename = filename.url();
823 
824  // Return
825  return;
826 }
827 
828 
829 /***********************************************************************//**
830  * @brief Save Pulse Height Analyzer spectrum
831  *
832  * @param[in] filename File name.
833  * @param[in] clobber Overwrite existing file?
834  *
835  * Saves the Pulse Height Analyzer spectrum and energy boundaries into a
836  * FITS file. If a file with the given @p filename does not yet exist it
837  * will be created. If the file exists it can be overwritten if the
838  * @p clobber flag is set to `true`. Otherwise an exception is thrown.
839  *
840  * The method will save two binary FITS tables into the FITS file: a
841  * `SPECTRUM` extension that contains the channel values of the Pulse
842  * Height Analyzer spectrum and an `EBOUNDS` extension that contains the
843  * energy boundaries for all channels.
844  ***************************************************************************/
845 void GPha::save(const GFilename& filename, const bool& clobber) const
846 {
847  // Create FITS file
848  GFits fits;
849 
850  // Write PHA into file
851  write(fits);
852 
853  // Save to file (without extension name since the requested extension
854  // may not yet exist in the file)
855  fits.saveto(filename.url(), clobber);
856 
857  // Store filename
858  m_filename = filename.url();
859 
860  // Return
861  return;
862 }
863 
864 
865 /***********************************************************************//**
866  * @brief Read Pulse Height Analyzer spectrum
867  *
868  * @param[in] fits File file.
869  *
870  * Reads the Pulse Height Analyzer spectrum from the `SPECTRUM` extension
871  * of the FITS file. If the file contains also an `EBOUNDS` extension the
872  * energy boundaries of all Pulse Height Analyzer channels are also loaded.
873  ***************************************************************************/
874 void GPha::read(const GFits& fits)
875 {
876  // Clear spectrum
877  clear();
878 
879  // Get PHA table
880  const GFitsTable& pha = *fits.table(gammalib::extname_pha);
881 
882  // Read PHA data
883  read(pha);
884 
885  // Optionally read energy boundaries
887 
888  // Get energy boundary table
890 
891  // Read energy boundaries
892  m_ebounds.read(ebounds);
893 
894  } // endif: had energy boundary table
895 
896  // Return
897  return;
898 }
899 
900 
901 /***********************************************************************//**
902  * @brief Read Pulse Height Analyzer spectrum
903  *
904  * @param[in] table FITS table.
905  *
906  * @exception GException::invalid_value
907  * Mismatch between PHA file and energy boundaries.
908  *
909  * Reads the Pulse Height Analyzer spectrum from a FITS table. The channel
910  * values are expected in the `COUNTS` column of the table. All other
911  * columns are ignored.
912  *
913  * See
914  * https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node5.html
915  * for details about the Pulse Height Analyzer spectrum format.
916  ***************************************************************************/
917 void GPha::read(const GFitsTable& table)
918 {
919  // Clear spectrum
920  clear();
921 
922  // Get data column
923  const GFitsTableCol* col_data = table["COUNTS"];
924  const GFitsTableCol* col_area = table["AREASCAL"];
925  const GFitsTableCol* col_back = table["BACKSCAL"];
926 
927  // Extract number of channels in FITS file
928  int length = col_data->nrows();
929 
930  // Check whether column length is consistent with energy boundaries
931  if (m_ebounds.size() > 0) {
932  if (m_ebounds.size() != length) {
933  std::string msg = "Mismatch between the "+gammalib::str(length)+
934  " channels in the PHA file and the "+
935  gammalib::str(m_ebounds.size())+" energy "
936  "boundaries. Please correct either the energy "
937  "boundaris or the PHA file.";
938  throw GException::invalid_value(G_READ, msg);
939  }
940  }
941 
942  // Initialize spectrum
943  alloc(length);
944 
945  // Copy data
946  for (int i = 0; i < length; ++i) {
947  m_counts[i] = col_data->real(i);
948  m_areascal[i] = col_area->real(i);
949  m_backscal[i] = col_back->real(i);
950  }
951 
952  // Read any additional columns
953  for (int icol = 0; icol < table.ncols(); ++icol) {
954 
955  // Fall through if the column is a standard column
956  std::string colname(table[icol]->name());
957  if ((colname == "CHANNEL") ||
958  (colname == "COUNTS") ||
959  (colname == "STAT_ERR") ||
960  (colname == "SYS_ERR") ||
961  (colname == "QUALITY") ||
962  (colname == "GROUPING") ||
963  (colname == "AREASCAL") ||
964  (colname == "BACKSCAL")) {
965  continue;
966  }
967 
968  // Get pointer to column
969  const GFitsTableCol* column = table[icol];
970 
971  // Set column vector
972  std::vector<double> coldata;
973  for (int i = 0; i < length; ++i) {
974  coldata.push_back(column->real(i));
975  }
976 
977  // Append column
978  append(colname, coldata);
979 
980  } // endfor: looped over all additional columns
981 
982  // Read file names
983  std::string backfile = (table.has_card("BACKFILE")) ? table.string("BACKFILE") : "";
984  std::string corrfile = (table.has_card("CORRFILE")) ? table.string("CORRFILE") : "";
985  std::string respfile = (table.has_card("RESPFILE")) ? table.string("RESPFILE") : "";
986  std::string ancrfile = (table.has_card("ANCRFILE")) ? table.string("ANCRFILE") : "";
987 
988  // Set file names
989  m_backfile = (gammalib::tolower(backfile) != "none") ? backfile : "";
990  m_corrfile = (gammalib::tolower(corrfile) != "none") ? corrfile : "";
991  m_respfile = (gammalib::tolower(respfile) != "none") ? respfile : "";
992  m_ancrfile = (gammalib::tolower(ancrfile) != "none") ? ancrfile : "";
993 
994  // Read energy band for observations
995  if (table.has_card("EMIN_OBS")) {
996  m_emin_obs.MeV(table.real("EMIN_OBS"));
997  }
998  if (table.has_card("EMAX_OBS")) {
999  m_emax_obs.MeV(table.real("EMAX_OBS"));
1000  }
1001 
1002  // Read keywords
1003  m_underflow = (table.has_card("UNDEFLOW")) ? table.real("UNDEFLOW") : 0.0;
1004  m_overflow = (table.has_card("OVERFLOW")) ? table.real("OVERFLOW") : 0.0;
1005  m_outflow = (table.has_card("OUTFLOW")) ? table.real("OUTFLOW") : 0.0;
1006  m_exposure = (table.has_card("EXPOSURE")) ? table.real("EXPOSURE") : 0.0;
1007 
1008  // Store FITS header
1009  m_header = table.header();
1010 
1011  // Return
1012  return;
1013 }
1014 
1015 
1016 /***********************************************************************//**
1017  * @brief Write Pulse Height Analyzer spectrum
1018  *
1019  * @param[in] fits FITS file.
1020  *
1021  * Writes the Pulse Height Analyzer spectrum into `SPECTRUM` and `EBOUNDS`
1022  * extensions of the FITS file. Extensions with these names will be removed
1023  * from the FITS file before writing.
1024  *
1025  * The columns `CHANNEL`, `COUNTS`, `STAT_ERR`, `SYS_ERR`, `QUALITY`,
1026  * `GROUPING`, `AREASCAL`, and `BACKSCAL` will be written into the `SPECTRUM`
1027  * extension, but only the `CHANNEL` and `COUNTS` columns will be filled with
1028  * values. Note that the channels start from 1 in the Pulse Height Analyzer
1029  * spectrum.
1030  *
1031  * See
1032  * https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node5.html
1033  * for details about the PHA file format.
1034  ***************************************************************************/
1035 void GPha::write(GFits& fits) const
1036 {
1037  // Remove extensions if they exist already
1038  if (fits.contains(gammalib::extname_ebounds)) {
1040  }
1041  if (fits.contains(gammalib::extname_pha)) {
1043  }
1044 
1045  // Set column length
1046  int length = size();
1047 
1048  // Continue only if there are bins
1049  if (length > 0) {
1050 
1051  // Create new binary table
1052  GFitsBinTable hdu;
1053 
1054  // Allocate floating point vector columns
1055  GFitsTableShortCol col_chan("CHANNEL", length);
1056  GFitsTableFloatCol col_data("COUNTS", length);
1057  GFitsTableFloatCol col_stat("STAT_ERR", length);
1058  GFitsTableFloatCol col_syst("SYS_ERR", length);
1059  GFitsTableShortCol col_qual("QUALITY", length);
1060  GFitsTableShortCol col_grpg("GROUPING", length);
1061  GFitsTableFloatCol col_area("AREASCAL", length);
1062  GFitsTableFloatCol col_back("BACKSCAL", length);
1063 
1064  // Fill columns
1065  for (int i = 0; i < length; ++i) {
1066  col_chan(i) = i+1; // Channels start at 1
1067  col_data(i) = float(m_counts[i]);
1068  col_stat(i) = float(std::sqrt(std::abs(m_counts[i])));
1069  col_grpg(i) = 1;
1070  col_area(i) = float(m_areascal[i]);
1071  col_back(i) = float(m_backscal[i]);
1072  }
1073 
1074  // Set table attributes
1076 
1077  // Append columns to table
1078  hdu.append(col_chan);
1079  hdu.append(col_data);
1080  hdu.append(col_stat);
1081  hdu.append(col_syst);
1082  hdu.append(col_qual);
1083  hdu.append(col_grpg);
1084  hdu.append(col_area);
1085  hdu.append(col_back);
1086 
1087  // Append any additional columns
1088  for (int icol = 0; icol < columns(); ++icol) {
1089 
1090  // Allocate floating point vector columns
1091  GFitsTableFloatCol column(m_colnames[icol], length);
1092 
1093  // Fill columns
1094  for (int i = 0; i < length; ++i) {
1095  column(i) = (float)m_coldata[icol][i];
1096  }
1097 
1098  // Append columns to table
1099  hdu.append(column);
1100 
1101  } // endfor: looped over all additional columns
1102 
1103  // Set file names
1104  std::string backfile = (m_backfile.empty()) ? "none" : m_backfile;
1105  std::string corrfile = (m_corrfile.empty()) ? "none" : m_corrfile;
1106  std::string respfile = (m_respfile.empty()) ? "none" : m_respfile;
1107  std::string ancrfile = (m_ancrfile.empty()) ? "none" : m_ancrfile;
1108 
1109  // Write header keywords
1110  hdu.card("TELESCOP", "unknown", "Telescope");
1111  hdu.card("INSTRUME", "unknown", "Instrument");
1112  hdu.card("FILTER", "none", "Filter");
1113  hdu.card("EXPOSURE", m_exposure, "[s] Deadtime corrected exposure time");
1114  hdu.card("BACKFILE", backfile, "Associated background filename");
1115  hdu.card("CORRFILE", corrfile, "Associated correction filename");
1116  hdu.card("CORRSCAL", 1.0, "Correction scaling factor");
1117  hdu.card("RESPFILE", respfile, "Associated RMF filename");
1118  hdu.card("ANCRFILE", ancrfile, "Associated ARF filename");
1119  hdu.card("HDUCLASS", "OGIP", "Format conforms to OGIP standard");
1120  hdu.card("HDUCLAS1", "SPECTRUM", "PHA dataset");
1121  hdu.card("HDUVERS", "1.2.1", "Version of the file format");
1122  hdu.card("POISSERR", true, "Poisson errors to be assumed");
1123  hdu.card("CHANTYPE", "PI", "Channel type");
1124  hdu.card("DETCHANS", size(), "Total number of possible channels");
1125  hdu.card("EMIN_OBS", m_emin_obs.MeV(), "[MeV] Minimum energy of observation");
1126  hdu.card("EMAX_OBS", m_emax_obs.MeV(), "[MeV] Maximum energy of observation");
1127  hdu.card("UNDEFLOW", m_underflow, "Number of underflowing events");
1128  hdu.card("OVERFLOW", m_overflow, "Number of overflowing events");
1129  hdu.card("OUTFLOW", m_outflow, "Number of outflowing events");
1130 
1131  // Write additional header keywords (do not overwrite any existing
1132  // keywords, except of "TELESCOP", "INSTRUME" and "FILTER")
1133  for (int i = 0; i < m_header.size(); ++i) {
1134  if ((!hdu.has_card(m_header[i].keyname())) ||
1135  (m_header[i].keyname() == "TELESCOP") ||
1136  (m_header[i].keyname() == "INSTRUME") ||
1137  (m_header[i].keyname() == "FILTER")) {
1138  hdu.card(m_header[i]);
1139  }
1140  }
1141 
1142  // Append HDU to FITS file
1143  fits.append(hdu);
1144 
1145  // Optionally append energy boundaries
1146  if (m_ebounds.size() > 0) {
1147  m_ebounds.write(fits);
1148  }
1149 
1150  } // endif: there were data to write
1151 
1152  // Return
1153  return;
1154 }
1155 
1156 
1157 /***********************************************************************//**
1158  * @brief Print Pulse Height Analyzer spectrum
1159  *
1160  * @param[in] chatter Chattiness.
1161  * @return String containing Pulse Height Analyzer spectrum information.
1162  ***************************************************************************/
1163 std::string GPha::print(const GChatter& chatter) const
1164 {
1165  // Initialise result string
1166  std::string result;
1167 
1168  // Continue only if chatter is not silent
1169  if (chatter != SILENT) {
1170 
1171  // Append header
1172  result.append("=== GPha ===");
1173 
1174  // Append energy boundary information
1175  result.append("\n"+gammalib::parformat("Exposure"));
1176  result.append(gammalib::str(exposure())+" s");
1177  result.append("\n"+gammalib::parformat("Number of bins"));
1178  result.append(gammalib::str(size()));
1179  result.append("\n"+gammalib::parformat("Energy range"));
1180  if (m_ebounds.size() > 0) {
1181  result.append(m_ebounds.emin().print());
1182  result.append(" - ");
1183  result.append(m_ebounds.emax().print());
1184  }
1185  else {
1186  result.append("not specified");
1187  }
1188 
1189  // Append observation energy range
1190  result.append("\n"+gammalib::parformat("Observation energy range"));
1191  if (m_emax_obs > m_emin_obs) {
1192  result.append(m_emin_obs.print());
1193  result.append(" - ");
1194  result.append(m_emax_obs.print());
1195  }
1196  else {
1197  result.append("not specified");
1198  }
1199 
1200  // Append content information
1201  result.append("\n"+gammalib::parformat("Total number of counts"));
1202  result.append(gammalib::str(counts()));
1203  result.append("\n"+gammalib::parformat("Underflow counts"));
1204  result.append(gammalib::str(m_underflow));
1205  result.append("\n"+gammalib::parformat("Overflow counts"));
1206  result.append(gammalib::str(m_overflow));
1207  result.append("\n"+gammalib::parformat("Outflow counts"));
1208  result.append(gammalib::str(m_outflow));
1209 
1210  } // endif: chatter was not silent
1211 
1212  // Return result
1213  return result;
1214 }
1215 
1216 
1217 /*==========================================================================
1218  = =
1219  = Private methods =
1220  = =
1221  ==========================================================================*/
1222 
1223 /***********************************************************************//**
1224  * @brief Initialise class members
1225  ***************************************************************************/
1227 {
1228  // Initialise members
1229  m_filename.clear();
1230  m_ebounds.clear();
1231  m_counts.clear();
1232  m_areascal.clear();
1233  m_backscal.clear();
1234  m_colnames.clear();
1235  m_coldata.clear();
1236  m_emin_obs.clear();
1237  m_emax_obs.clear();
1238  m_underflow = 0.0;
1239  m_overflow = 0.0;
1240  m_outflow = 0.0;
1241  m_exposure = 0.0;
1242  m_backfile.clear();
1243  m_corrfile.clear();
1244  m_respfile.clear();
1245  m_ancrfile.clear();
1246  m_header.clear();
1247 
1248  // Return
1249  return;
1250 }
1251 
1252 
1253 /***********************************************************************//**
1254  * @brief Copy class members
1255  *
1256  * @param[in] pha Pulse Height Analyzer spectrum.
1257  ***************************************************************************/
1258 void GPha::copy_members(const GPha& pha)
1259 {
1260  // Copy members
1261  m_filename = pha.m_filename;
1262  m_counts = pha.m_counts;
1263  m_areascal = pha.m_areascal;
1264  m_backscal = pha.m_backscal;
1265  m_colnames = pha.m_colnames;
1266  m_coldata = pha.m_coldata;
1267  m_emin_obs = pha.m_emin_obs;
1268  m_emax_obs = pha.m_emax_obs;
1269  m_underflow = pha.m_underflow;
1270  m_overflow = pha.m_overflow;
1271  m_outflow = pha.m_outflow;
1272  m_ebounds = pha.m_ebounds;
1273  m_exposure = pha.m_exposure;
1274  m_backfile = pha.m_backfile;
1275  m_corrfile = pha.m_corrfile;
1276  m_respfile = pha.m_respfile;
1277  m_ancrfile = pha.m_ancrfile;
1278  m_header = pha.m_header;
1279 
1280  // Return
1281  return;
1282 }
1283 
1284 
1285 /***********************************************************************//**
1286  * @brief Delete class members
1287  ***************************************************************************/
1289 {
1290  // Return
1291  return;
1292 }
1293 
1294 
1295 /***********************************************************************//**
1296  * @brief Allocate spectrum
1297  *
1298  * @param[in] size Size of spectrum.
1299  *
1300  * Allocates memory for a spectrum of given @p size and sets number of counts
1301  * to zero and area and background scaling factors to one.
1302  ***************************************************************************/
1303 void GPha::alloc(const int& size)
1304 {
1305  // Initialise attributes
1306  m_underflow = 0.0;
1307  m_overflow = 0.0;
1308  m_outflow = 0.0;
1309  m_exposure = 0.0;
1310 
1311  // Initialize vectors
1312  m_counts.assign(size, 0.0);
1313  m_areascal.assign(size, 1.0);
1314  m_backscal.assign(size, 1.0);
1315 
1316  // Return
1317  return;
1318 }
1319 
1320 
1321 /***********************************************************************//**
1322  * @brief Returns index of additional vector column
1323  *
1324  * @param[in] colname Vector column name.
1325  * @return Vector column index.
1326  *
1327  * Returns index of additional vector column. Returns -1 if the vector column
1328  * has not found.
1329  ***************************************************************************/
1330 int GPha::column_index(const std::string& colname) const
1331 {
1332  // Initialise index
1333  int index(-1);
1334 
1335  // Search vector column name
1336  for (int i = 0; i < columns(); ++i) {
1337  if (m_colnames[i] == colname) {
1338  index = i;
1339  break;
1340  }
1341  }
1342 
1343  // Return index
1344  return index;
1345 }
const int & ncols(void) const
Return number of columns in table.
Definition: GFitsTable.hpp:134
const double & exposure(void) const
Return exposure time.
Definition: GPha.hpp:355
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
const std::string & ancrfile(void) const
Return Ancilliary Response File name.
Definition: GPha.hpp:529
bool contains(const int &extno) const
Check if HDU exists in FITS file.
Definition: GFits.hpp:282
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
GPha * clone(void) const
Clone object.
Definition: GPha.cpp:465
void append(const std::string &name, const std::vector< double > &column)
Append additional column to spectrum.
Definition: GPha.cpp:574
virtual ~GPha(void)
Destructor.
Definition: GPha.cpp:159
Energy value class definition.
double m_overflow
Number of overflowing events.
Definition: GPha.hpp:164
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
double m_exposure
Deadtime corr. exp. time (sec)
Definition: GPha.hpp:166
const GEbounds & ebounds(void) const
Return energy boundaries.
Definition: GPha.hpp:282
double m_underflow
Number of underflowing events.
Definition: GPha.hpp:163
double counts(void) const
Number of counts in spectrum.
Definition: GPha.cpp:703
#define G_READ
Definition: GPha.cpp:53
void nrows(const int &nrows)
Set number of rows in column.
#define G_OPERATOR
Definition: GPha.cpp:45
GNdarray counts_spectrum(void) const
Get number of counts in spectrum as GNdarray.
Definition: GPha.cpp:724
GFilename m_filename
Filename of origin.
Definition: GPha.hpp:155
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition: GFitsTable.hpp:147
void backscal(const int &index, const double &backscal)
Set background scaling factor.
Definition: GPha.cpp:658
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
std::vector< double > m_areascal
Area scaling.
Definition: GPha.hpp:157
void read(const GFitsTable &table)
Read energy boundaries from FITS table.
Definition: GEbounds.cpp:776
#define G_AT2
Definition: GPha.cpp:47
Gammalib tools definition.
FITS table float column class interface definition.
GPha & operator*=(const double &scale)
Scale spectrum values.
Definition: GPha.cpp:343
FITS file class.
Definition: GFits.hpp:63
const GFitsHeader & header(void) const
Return extension header.
Definition: GFitsHDU.hpp:205
FITS file class interface definition.
void fill(const GEnergy &energy, const double &value=1.0)
Fill spectrum with a value.
Definition: GPha.cpp:771
void init_members(void)
Initialise class members.
Definition: GPha.cpp:1226
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:831
void load(const GFilename &filename)
Load Pulse Height Analyzer spectrum.
Definition: GPha.cpp:806
int index(const GEnergy &eng) const
Returns energy bin index for a given energy.
Definition: GEbounds.cpp:963
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1267
double m_outflow
Number of outflowing events.
Definition: GPha.hpp:165
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
#define G_BACKSCAL_SET
Definition: GPha.cpp:51
#define G_OPERATOR_PLUS
Definition: GPha.cpp:43
GEbounds m_ebounds
Energy boundaries.
Definition: GPha.hpp:167
XSPEC Pulse Height Analyzer class definition.
int column_index(const std::string &colname) const
Returns index of additional vector column.
Definition: GPha.cpp:1330
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
#define G_OPERATOR_MINUS
Definition: GPha.cpp:44
std::vector< std::vector< double > > m_coldata
Additional column data.
Definition: GPha.hpp:160
Energy boundaries container class.
Definition: GEbounds.hpp:60
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:848
const std::string & corrfile(void) const
Return correction file name.
Definition: GPha.hpp:471
std::string m_respfile
RMF file.
Definition: GPha.hpp:170
const std::string & respfile(void) const
Return Redistribution Matrix File name.
Definition: GPha.hpp:500
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
Filename class.
Definition: GFilename.hpp:62
N-dimensional array class interface definition.
std::string m_corrfile
Correction file.
Definition: GPha.hpp:169
GPha & operator=(const GPha &pha)
Assignment operator.
Definition: GPha.cpp:181
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
Abstract interface for FITS table column.
double & operator[](const int &index)
Return content of spectral bin.
Definition: GPha.hpp:196
void free_members(void)
Delete class members.
Definition: GPha.cpp:1288
const std::string extname_ebounds
Definition: GEbounds.hpp:44
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
int size(void) const
Return number of cards in header.
GChatter
Definition: GTypemaps.hpp:33
#define G_APPEND
Definition: GPha.cpp:48
N-dimensional array class.
Definition: GNdarray.hpp:44
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
const std::string extname_pha
Definition: GPha.hpp:45
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:269
void clear(void)
Clear object.
Definition: GPha.cpp:447
GPha & operator/=(const double &scale)
Divide spectrum values.
Definition: GPha.cpp:369
const std::string & backfile(void) const
Return background file name.
Definition: GPha.hpp:442
GPha & operator+=(const GPha &pha)
Add spectrum.
Definition: GPha.cpp:231
void write(GFits &fits) const
Write Pulse Height Analyzer spectrum.
Definition: GPha.cpp:1035
void save(const GFilename &filename, const bool &clobber=false) const
Save Pulse Height Analyzer spectrum.
Definition: GPha.cpp:845
std::vector< double > m_counts
Counts data.
Definition: GPha.hpp:156
std::string m_backfile
Background file.
Definition: GPha.hpp:168
GEnergy m_emin_obs
Minimum energy of observation.
Definition: GPha.hpp:161
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
#define G_AREASCAL_GET
Definition: GPha.cpp:50
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
void copy_members(const GPha &pha)
Copy class members.
Definition: GPha.cpp:1258
virtual double real(const int &row, const int &inx=0) const =0
std::string print(const GChatter &chatter=NORMAL) const
Print Pulse Height Analyzer spectrum.
Definition: GPha.cpp:1163
int size(void) const
Return number of bins in spectrum.
Definition: GPha.hpp:254
#define G_AT1
Definition: GPha.cpp:46
FITS binary table class.
void read(const GFits &fits)
Read Pulse Height Analyzer spectrum.
Definition: GPha.cpp:874
GPha(void)
Void constructor.
Definition: GPha.cpp:71
std::string string(const std::string &keyname) const
Return card value as string.
Definition: GFitsHDU.hpp:410
Exception handler interface definition.
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition: GFits.cpp:665
GFitsHeader m_header
FITS header cards.
Definition: GPha.hpp:172
FITS binary table class definition.
std::string tolower(const std::string &s)
Convert string to lower case.
Definition: GTools.cpp:834
FITS table short integer column class interface definition.
double & at(const int &index)
Return content of spectral bin.
Definition: GPha.cpp:482
FITS table float column.
void areascal(const int &index, const double &areascal)
Set area scaling factor.
Definition: GPha.cpp:607
void alloc(const int &size)
Allocate spectrum.
Definition: GPha.cpp:1303
FITS table short integer column.
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
std::string m_ancrfile
ARF file.
Definition: GPha.hpp:171
void clear(void)
Clear header.
std::vector< std::string > m_colnames
Additional column names.
Definition: GPha.hpp:159
std::vector< double > m_backscal
Background scaling.
Definition: GPha.hpp:158
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
int columns(void) const
Return number of additional columns.
Definition: GPha.hpp:268
Pulse Height Analyzer class.
Definition: GPha.hpp:64
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
GNdarray backscal_spectrum(void) const
Get background scaling factors as GNdarray.
Definition: GPha.cpp:747
GPha & operator-=(const GPha &pha)
Subtract spectrum.
Definition: GPha.cpp:297
#define G_AREASCAL_SET
Definition: GPha.cpp:49
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
GEnergy m_emax_obs
Minimum energy of observation.
Definition: GPha.hpp:162
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
#define G_BACKSCAL_GET
Definition: GPha.cpp:52
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
FITS table abstract base class interface definition.