GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GArf.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GArf.cpp - XSPEC Auxiliary Response File class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2013-2021 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 GArf.cpp
23  * @brief XSPEC Auxiliary Response File class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GArf.hpp"
32 #include "GException.hpp"
33 #include "GFilename.hpp"
34 #include "GTools.hpp"
35 #include "GFits.hpp"
36 #include "GFitsTable.hpp"
37 #include "GFitsBinTable.hpp"
38 #include "GFitsTableShortCol.hpp"
39 #include "GFitsTableFloatCol.hpp"
40 
41 /* __ Method name definitions ____________________________________________ */
42 #define G_OPERATOR_PLUS "GArf::operator+=(GArf&)"
43 #define G_OPERATOR_MINUS "GArf::operator-=(GArf&)"
44 #define G_OPERATOR "GArf::operator[](std::string&)"
45 #define G_OPERATOR2 "GArf::operator()(std::string&, GEnergy&)"
46 #define G_AT1 "GArf::at(int&)"
47 #define G_AT2 "GArf::at(int&, int&)"
48 #define G_APPEND "GArf::append(std::string&, std::vector<double>&)"
49 
50 /* __ Macros _____________________________________________________________ */
51 
52 /* __ Coding definitions _________________________________________________ */
53 
54 /* __ Debug definitions __________________________________________________ */
55 
56 
57 /*==========================================================================
58  = =
59  = Constructors/destructors =
60  = =
61  ==========================================================================*/
62 
63 /***********************************************************************//**
64  * @brief Void constructor
65  ***************************************************************************/
67 {
68  // Initialise members
69  init_members();
70 
71  // Return
72  return;
73 }
74 
75 
76 /***********************************************************************//**
77  * @brief File constructor
78  *
79  * @param[in] filename File name.
80  ***************************************************************************/
81 GArf::GArf(const GFilename& filename)
82 {
83  // Initialise members
84  init_members();
85 
86  // Load ARF file
87  load(filename);
88 
89  // Return
90  return;
91 }
92 
93 
94 /***********************************************************************//**
95  * @brief Energy boundary constructor
96  *
97  * @param[in] ebds Energy boundaries.
98  ***************************************************************************/
99 GArf::GArf(const GEbounds& ebds)
100 {
101  // Initialise members
102  init_members();
103 
104  // Set energy boundaries
105  m_ebounds = ebds;
106 
107  // Set log true energy node array
108  set_logetrue();
109 
110  // Initialize spectral response
111  m_specresp.assign(ebds.size(), 0.0);
112 
113  // Return
114  return;
115 }
116 
117 
118 /***********************************************************************//**
119  * @brief Copy constructor
120  *
121  * @param[in] arf Auxiliary Response File.
122  ***************************************************************************/
123 GArf::GArf(const GArf& arf)
124 {
125  // Initialise members
126  init_members();
127 
128  // Copy members
129  copy_members(arf);
130 
131  // Return
132  return;
133 }
134 
135 
136 /***********************************************************************//**
137  * @brief Destructor
138  ***************************************************************************/
140 {
141  // Free members
142  free_members();
143 
144  // Return
145  return;
146 }
147 
148 
149 /*==========================================================================
150  = =
151  = Operators =
152  = =
153  ==========================================================================*/
154 
155 /***********************************************************************//**
156  * @brief Assignment operator
157  *
158  * @param[in] arf Auxiliary Response File.
159  * @return Auxiliary Response File.
160  ***************************************************************************/
162 {
163  // Execute only if object is not identical
164  if (this != &arf) {
165 
166  // Free members
167  free_members();
168 
169  // Initialise members
170  init_members();
171 
172  // Copy members
173  copy_members(arf);
174 
175  } // endif: object was not identical
176 
177  // Return
178  return *this;
179 }
180 
181 
182 /***********************************************************************//**
183  * @brief Add Auxiliary Response File
184  *
185  * @param[in] arf Auxiliary Response File.
186  * @return Sum of Auxiliary Response Files.
187  *
188  * @exception GException::invalid_value
189  * Incompatible Auxiliary Response Files.
190  *
191  * Adds the ARF values of an Auxiliary Response File to the current values.
192  *
193  * The operator only works if the provide Auxiliary Response File has the
194  * same energy binning than the current Auxiliary Response File.
195  ***************************************************************************/
197 {
198  // Throw an exception if the ARF are not compatible
199  if (this->ebounds() != arf.ebounds()) {
200  std::string msg = "Incompatible energy binning of Auxiliary "
201  "Response File.";
203  }
204 
205  // Add ARF values
206  for (int i = 0; i < this->size(); ++i) {
207  m_specresp[i] += arf.m_specresp[i];
208  }
209 
210  // Return
211  return *this;
212 }
213 
214 
215 /***********************************************************************//**
216  * @brief Subtract Auxiliary Response File
217  *
218  * @param[in] arf Auxiliary Response File.
219  * @return Difference of Auxiliary Response File.
220  *
221  * @exception GException::invalid_value
222  * Incompatible Auxiliary Response Files.
223  *
224  * Subtracts the ARF values of an Auxiliary Response File from the current
225  * values.
226  *
227  * The operator only works if the provide Auxiliary Response File has the
228  * same energy binning than the current Auxiliary Response File.
229  ***************************************************************************/
231 {
232  // Throw an exception if the ARF are not compatible
233  if (this->ebounds() != arf.ebounds()) {
234  std::string msg = "Incompatible energy binning of Auxiliary "
235  "Response File.";
237  }
238 
239  // Subtract ARF values
240  for (int i = 0; i < this->size(); ++i) {
241  m_specresp[i] -= arf.m_specresp[i];
242  }
243 
244  // Return
245  return *this;
246 }
247 
248 
249 /***********************************************************************//**
250  * @brief Scale Auxiliary Response File values
251  *
252  * @param[in] scale Scale factor.
253  * @return Scaled Auxiliary Response File.
254  *
255  * Multiplies the values of the Auxiliary Response File with a scale factor.
256  ***************************************************************************/
257 GArf& GArf::operator*=(const double& scale)
258 {
259  // Scale ARF values
260  for (int i = 0; i < this->size(); ++i) {
261  m_specresp[i] *= scale;
262  }
263 
264  // Return
265  return *this;
266 }
267 
268 
269 /***********************************************************************//**
270  * @brief Divide Auxiliary Response File values
271  *
272  * @param[in] scale Division factor.
273  * @return Divided Auxiliary Response File.
274  *
275  * Divides the values of the Auxiliary Response File by a division factor.
276  ***************************************************************************/
277 GArf& GArf::operator/=(const double& scale)
278 {
279  // Divide ARF values
280  for (int i = 0; i < this->size(); ++i) {
281  m_specresp[i] /= scale;
282  }
283 
284  // Return
285  return *this;
286 }
287 
288 
289 /***********************************************************************//**
290  * @brief Return additional vector column
291  *
292  * @param[in] colname Vector column name.
293  * @return Vector column.
294  *
295  * Returns reference to additional vector column.
296  ***************************************************************************/
297 std::vector<double>& GArf::operator[](const std::string& colname)
298 {
299  // Determine index of additional column (-1 if not found)
300  int index = column_index(colname);
301 
302  // Throw exception if index not found
303  if (index == -1) {
304  std::string msg = "Could not find additional column with name \""+
305  colname+"\".";
307  }
308 
309  // Return reference to vector column
310  return (m_coldata[index]);
311 }
312 
313 
314 /***********************************************************************//**
315  * @brief Return additional vector column (const version)
316  *
317  * @param[in] colname Vector column name.
318  * @return Vector column.
319  *
320  * Returns reference to additional vector column.
321  ***************************************************************************/
322 const std::vector<double>& GArf::operator[](const std::string& colname) const
323 {
324  // Determine index of additional column (-1 if not found)
325  int index = column_index(colname);
326 
327  // Throw exception if index not found
328  if (index == -1) {
329  std::string msg = "Could not find additional column with name \""+
330  colname+"\".";
332  }
333 
334  // Return reference to vector column
335  return (m_coldata[index]);
336 }
337 
338 
339 /***********************************************************************//**
340  * @brief Return vector column content as function of energy
341  *
342  * @param[in] colname Vector column name.
343  * @param[in] energy Energy.
344  * @return Vector column content.
345  *
346  * Returns content of additional vector column as function of @p energy.
347  ***************************************************************************/
348 double GArf::operator()(const std::string& colname,
349  const GEnergy& energy) const
350 {
351  // Determine index of additional column (-1 if not found)
352  int index = column_index(colname);
353 
354  // Throw exception if index not found
355  if (index == -1) {
356  std::string msg = "Could not find additional column with name \""+
357  colname+"\".";
359  }
360 
361  // Get reference to vector column
362  const std::vector<double>& column = m_coldata[index];
363 
364  // Perform log-log interpolation of vector column content
365  m_logetrue.set_value(energy.log10TeV());
366  double wgt_left = m_logetrue.wgt_left();
367  double wgt_right = m_logetrue.wgt_right();
368  double value_left = column[m_logetrue.inx_left()];
369  double value_right = column[m_logetrue.inx_right()];
370  double value = 0.0;
371  if (value_left > 0.0 && value_right > 0.0) {
372  value = std::exp(wgt_left * std::log(value_left) +
373  wgt_right * std::log(value_right));
374  }
375  if (value < 0.0) {
376  value = 0.0;
377  }
378 
379  // Return value
380  return value;
381 }
382 
383 
384 /*==========================================================================
385  = =
386  = Public methods =
387  = =
388  ==========================================================================*/
389 
390 /***********************************************************************//**
391  * @brief Clear object.
392  *
393  * Reset object to a clean initial state.
394  ***************************************************************************/
395 void GArf::clear(void)
396 {
397  // Free memory and initialise members
398  free_members();
399  init_members();
400 
401  // Return
402  return;
403 }
404 
405 
406 /***********************************************************************//**
407  * @brief Clone object
408  *
409  * @return Pointer to Auxiliary Response File.
410  ***************************************************************************/
411 GArf* GArf::clone(void) const
412 {
413  // Clone client
414  return new GArf(*this);
415 }
416 
417 
418 /***********************************************************************//**
419  * @brief Return content of spectral bin
420  *
421  * @param[in] index Spectral bin index [0,...,size()[.
422  *
423  * @exception GException::out_of_range
424  * Spectral bin index is out of range.
425  *
426  * Returns reference to content of spectral bin with specified @p index.
427  ***************************************************************************/
428 double& GArf::at(const int& index)
429 {
430  // Throw an exception if index is out of range
431  if (index < 0 || index >= size()) {
432  throw GException::out_of_range(G_AT1, "Spectral bin index",
433  index, size());
434  }
435 
436  // Return reference
437  return (m_specresp[index]);
438 }
439 
440 
441 /***********************************************************************//**
442  * @brief Return content of spectral bin (const version)
443  *
444  * @param[in] index Spectral bin index [0,...,size()[.
445  *
446  * @exception GException::out_of_range
447  * Spectral bin index is out of range.
448  * *
449  * Returns reference to content of spectral bin with specified @p index.
450  ***************************************************************************/
451 const double& GArf::at(const int& index) const
452 {
453  // Throw an exception if index is out of range
454  if (index < 0 || index >= size()) {
455  throw GException::out_of_range(G_AT1, "Spectral bin index",
456  index, size());
457  }
458 
459  // Return reference
460  return (m_specresp[index]);
461 }
462 
463 
464 /***********************************************************************//**
465  * @brief Return content of additional columns
466  *
467  * @param[in] index Bin index [0,...,size()-1].
468  * @param[in] col Columns index [0,...,columns()-1].
469  *
470  * @exception GException::out_of_range
471  * Spectral bin or column index is out of range.
472  *
473  * Returns reference to content of additional columns.
474  ***************************************************************************/
475 double& GArf::at(const int& index, const int& col)
476 {
477  // Throw an exception if bin or column index is out of range
478  if (index < 0 || index >= size()) {
479  throw GException::out_of_range(G_AT2, "Bin index", index, size());
480  }
481  if (col < 0 || col >= columns()) {
482  throw GException::out_of_range(G_AT2, "Column index", col, columns());
483  }
484 
485  // Return reference
486  return (m_coldata[col][index]);
487 }
488 
489 
490 /***********************************************************************//**
491  * @brief Return content of additional columns (const version)
492  *
493  * @param[in] index Bin index [0,...,size()-1].
494  * @param[in] col Columns index [0,...,columns()-1].
495  *
496  * @exception GException::out_of_range
497  * Spectral bin or column index is out of range.
498  *
499  * Returns reference to content of additional columns.
500  ***************************************************************************/
501 const double& GArf::at(const int& index, const int& col) const
502 {
503  // Throw an exception if bin or column index is out of range
504  if (index < 0 || index >= size()) {
505  throw GException::out_of_range(G_AT2, "Bin index", index, size());
506  }
507  if (col < 0 || col >= columns()) {
508  throw GException::out_of_range(G_AT2, "Column index", col, columns());
509  }
510 
511  // Return reference
512  return (m_coldata[col][index]);
513 }
514 
515 
516 /***********************************************************************//**
517  * @brief Append additional column to spectrum
518  *
519  * @param[in] name Additional column name.
520  * @param[in] column Additional column data.
521  ***************************************************************************/
522 void GArf::append(const std::string& name, const std::vector<double>& column)
523 {
524  // Throw an exception if the number of elements in the column does not
525  // correspond to the size of the spectrum
526  if (column.size() != size()) {
527  std::string msg = "Size of column "+gammalib::str(column.size())+
528  " is incompatible with size of spectrum "+
529  gammalib::str(size())+".";
531  }
532 
533  // Append column name
534  m_colnames.push_back(name);
535 
536  // Append column data
537  m_coldata.push_back(column);
538 
539  // Return
540  return;
541 }
542 
543 
544 /***********************************************************************//**
545  * @brief Load Auxiliary Response File
546  *
547  * @param[in] filename File name.
548  *
549  * Loads the Auxiliary Response File from the `SPECRESP` extension of the
550  * FITS file.
551  ***************************************************************************/
552 void GArf::load(const GFilename& filename)
553 {
554  // Clear response
555  clear();
556 
557  // Open FITS file (without extension name as the user is not allowed
558  // to modify the extension names)
559  GFits fits(filename.url());
560 
561  // Read ARF data
562  read(fits);
563 
564  // Close FITS file
565  fits.close();
566 
567  // Store filename
568  m_filename = filename.url();
569 
570  // Return
571  return;
572 }
573 
574 
575 /***********************************************************************//**
576  * @brief Save Auxiliary Response File
577  *
578  * @param[in] filename File name.
579  * @param[in] clobber Overwrite existing file?
580  *
581  * Saves the Auxiliary Response File into a FITS file. If a file with the
582  * given @p filename does not yet exist it will be created. If the file
583  * exists it can be overwritten if the @p clobber flag is set to `true`.
584  * Otherwise an exception is thrown.
585  *
586  * The method will save the `SPECRESP` binary FITS table into the FITS file
587  * that contains the values of the Auxiliary Response File.
588  ***************************************************************************/
589 void GArf::save(const GFilename& filename, const bool& clobber) const
590 {
591  // Create FITS file
592  GFits fits;
593 
594  // Write ARF into file
595  write(fits);
596 
597  // Save to file (without extension name since the requested extension
598  // may not yet exist in the file)
599  fits.saveto(filename.url(), clobber);
600 
601  // Store filename
602  m_filename = filename.url();
603 
604  // Return
605  return;
606 }
607 
608 
609 /***********************************************************************//**
610  * @brief Read Auxiliary Response File
611  *
612  * @param[in] fits FITS file.
613  *
614  * Loads the Auxiliary Response File from the `SPECRESP` extension of the
615  * FITS file.
616  ***************************************************************************/
617 void GArf::read(const GFits& fits)
618 {
619  // Clear response
620  clear();
621 
622  // Get ARF table
623  const GFitsTable& table = *fits.table(gammalib::extname_arf);
624 
625  // Read ARF data
626  read(table);
627 
628  // Return
629  return;
630 }
631 
632 
633 /***********************************************************************//**
634  * @brief Read Auxiliary Response File
635  *
636  * @param[in] table ARF FITS table.
637  *
638  * Reads the Auxiliary Response File from a FITS table. The true energy
639  * boundaries are expected in the `ENERG_LO` and `ENERG_HI` columns, the
640  * response information is expected in the `SPECRESP` column.
641  *
642  * The method will analyze the unit of the `SPECRESP` column, and if either
643  * `m**2`, `m^2` or `m2` are encountered, multiply the values of the column
644  * by \f$10^4\f$ to convert the response into units of \f$cm^2\f$. Units of
645  * the `ENERG_LO` and `ENERG_HI` columns are also interpreted for conversion.
646  *
647  * See
648  * http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html#tth_sEc4
649  * for details about the Auxiliary Response File format.
650  ***************************************************************************/
651 void GArf::read(const GFitsTable& table)
652 {
653  // Clear members
654  clear();
655 
656  // Get pointer to data columns
657  const GFitsTableCol* energy_lo = table["ENERG_LO"];
658  const GFitsTableCol* energy_hi = table["ENERG_HI"];
659  const GFitsTableCol* specresp = table["SPECRESP"];
660 
661  // Determine effective area conversion factor. Internal
662  // units are cm^2
663  std::string u_specresp =
665  double c_specresp = 1.0;
666  if (u_specresp == "m**2" || u_specresp == "m^2" || u_specresp == "m2") {
667  c_specresp = 10000.0;
668  }
669 
670  // Extract number of energy bins
671  int num = energy_lo->nrows();
672 
673  // Set energy bins
674  for (int i = 0; i < num; ++i) {
675 
676  // Append energy bin
677  GEnergy emin(energy_lo->real(i), energy_lo->unit());
678  GEnergy emax(energy_hi->real(i), energy_hi->unit());
679  m_ebounds.append(emin, emax);
680 
681  // Append effective area value
682  double aeff = specresp->real(i) * c_specresp;
683  m_specresp.push_back(aeff);
684 
685  } // endfor: looped over energy bins
686 
687  // Read any additional columns
688  for (int icol = 0; icol < table.ncols(); ++icol) {
689 
690  // Fall through if the column is a standard column
691  std::string colname(table[icol]->name());
692  if ((colname == "ENERG_LO") ||
693  (colname == "ENERG_HI") ||
694  (colname == "SPECRESP")) {
695  continue;
696  }
697 
698  // Get pointer to column
699  const GFitsTableCol* column = table[icol];
700 
701  // Set column vector
702  std::vector<double> coldata;
703  for (int i = 0; i < num; ++i) {
704  coldata.push_back(column->real(i));
705  }
706 
707  // Append column
708  append(colname, coldata);
709 
710  } // endfor: looped over all additional columns
711 
712  // Set log true energy node array
713  set_logetrue();
714 
715  // Store FITS header
716  m_header = table.header();
717 
718  // Return
719  return;
720 }
721 
722 
723 /***********************************************************************//**
724  * @brief Write Auxiliary Response File
725  *
726  * @param[in] fits FITS file.
727  *
728  * Writes the Auxiliary Response File into the `SPECRESP` extension of the
729  * FITS file. An existing extension with the same name will be removed from
730  * the FITS file before writing.
731  *
732  * The method writes the boundaries of the true energy bins into the
733  * `ENERG_LO` and `ENERG_HI` columns, response information will be written
734  * into the `SPECRESP` column.
735  *
736  * See
737  * http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html#tth_sEc4
738  * for details about the Auxiliary Response File format.
739  ***************************************************************************/
740 void GArf::write(GFits& fits) const
741 {
742  // If the FITS object contains already an extension with the same
743  // name then remove now this extension
744  if (fits.contains(gammalib::extname_arf)) {
746  }
747 
748  // Set column length
749  int length = size();
750 
751  // Continue only if there are bins
752  if (length > 0) {
753 
754  // Create new binary table
755  GFitsBinTable hdu;
756 
757  // Allocate floating point vector columns
758  GFitsTableFloatCol energy_lo("ENERG_LO", length);
759  GFitsTableFloatCol energy_hi("ENERG_HI", length);
760  GFitsTableFloatCol specresp("SPECRESP", length);
761 
762  // Fill columns
763  for (int i = 0; i < length; ++i) {
764  energy_lo(i) = (float)m_ebounds.emin(i).keV();
765  energy_hi(i) = (float)m_ebounds.emax(i).keV();
766  specresp(i) = (float)m_specresp[i];
767  }
768 
769  // Set column units
770  energy_lo.unit("keV");
771  energy_hi.unit("keV");
772  specresp.unit("cm**2");
773 
774  // Set table attributes
776 
777  // Append columns to table
778  hdu.append(energy_lo);
779  hdu.append(energy_hi);
780  hdu.append(specresp);
781 
782  // Append any additional columns
783  for (int icol = 0; icol < columns(); ++icol) {
784 
785  // Allocate floating point vector columns
786  GFitsTableFloatCol column(m_colnames[icol], length);
787 
788  // Fill columns
789  for (int i = 0; i < length; ++i) {
790  column(i) = (float)m_coldata[icol][i];
791  }
792 
793  // Append columns to table
794  hdu.append(column);
795 
796  } // endfor: looped over all additional columns
797 
798  // Write mandatory header keywords
799  hdu.card("TELESCOP", "unknown", "Telescope");
800  hdu.card("INSTRUME", "unknown", "Instrument");
801  hdu.card("FILTER", "none", "Filter");
802  hdu.card("HDUCLASS", "OGIP", "Format conforms to OGIP standard");
803  hdu.card("HDUCLAS1", "RESPONSE", "Extension contains response data");
804  hdu.card("HDUCLAS2", "SPECRESP", "Extension contains an ARF");
805  hdu.card("HDUVERS", "1.1.0", "Version of the file format");
806 
807  // Write additional header keywords (do not overwrite any existing
808  // keywords, except of "TELESCOP", "INSTRUME" and "FILTER")
809  for (int i = 0; i < m_header.size(); ++i) {
810  if ((!hdu.contains(m_header[i].keyname())) ||
811  (m_header[i].keyname() == "TELESCOP") ||
812  (m_header[i].keyname() == "INSTRUME") ||
813  (m_header[i].keyname() == "FILTER")) {
814  hdu.card(m_header[i]);
815  }
816  }
817 
818  // Append HDU to FITS file
819  fits.append(hdu);
820 
821  } // endif: there were data to write
822 
823  // Return
824  return;
825 }
826 
827 
828 /***********************************************************************//**
829  * @brief Print Auxiliary Response File
830  *
831  * @param[in] chatter Chattiness.
832  * @return String containing Auxiliary Response File information.
833  ***************************************************************************/
834 std::string GArf::print(const GChatter& chatter) const
835 {
836  // Initialise result string
837  std::string result;
838 
839  // Continue only if chatter is not silent
840  if (chatter != SILENT) {
841 
842  // Append header
843  result.append("=== GArf ===");
844 
845  // Append energy boundary information
846  result.append("\n"+gammalib::parformat("Number of bins"));
847  result.append(gammalib::str(m_ebounds.size()));
848  result.append("\n"+gammalib::parformat("Energy range"));
849  result.append(m_ebounds.emin().print());
850  result.append(" - ");
851  result.append(m_ebounds.emax().print());
852 
853  } // endif: chatter was not silent
854 
855  // Return result
856  return result;
857 }
858 
859 
860 /*==========================================================================
861  = =
862  = Private methods =
863  = =
864  ==========================================================================*/
865 
866 /***********************************************************************//**
867  * @brief Initialise class members
868  ***************************************************************************/
870 {
871  // Initialise members
872  m_filename.clear();
873  m_ebounds.clear();
874  m_logetrue.clear();
875  m_specresp.clear();
876  m_colnames.clear();
877  m_coldata.clear();
878  m_header.clear();
879 
880  // Return
881  return;
882 }
883 
884 
885 /***********************************************************************//**
886  * @brief Copy class members
887  *
888  * @param[in] arf Auxiliary Response File.
889  ***************************************************************************/
890 void GArf::copy_members(const GArf& arf)
891 {
892  // Copy members
893  m_filename = arf.m_filename;
894  m_ebounds = arf.m_ebounds;
895  m_logetrue = arf.m_logetrue;
896  m_specresp = arf.m_specresp;
897  m_colnames = arf.m_colnames;
898  m_coldata = arf.m_coldata;
899  m_header = arf.m_header;
900 
901  // Return
902  return;
903 }
904 
905 
906 /***********************************************************************//**
907  * @brief Delete class members
908  ***************************************************************************/
910 {
911  // Return
912  return;
913 }
914 
915 
916 /***********************************************************************//**
917  * @brief Set true energy node array
918  ***************************************************************************/
920 {
921  // Clear node array
922  m_logetrue.clear();
923 
924  // Continue only if there are true energies in Arf
925  int netrue = m_ebounds.size();
926  if (netrue > 0) {
927 
928  // Reserve space in node array
929  m_logetrue.reserve(netrue);
930 
931  // Append all log mean energies to node array
932  for (int i = 0; i < netrue; ++i) {
933 
934  // Get log mean of true energy in TeV
935  double logE = m_ebounds.elogmean(i).log10TeV();
936 
937  // Append energy to node array
938  m_logetrue.append(logE);
939 
940  } // endfor: appended log mean energies
941 
942  } // endif: there were true energies in Arf
943 
944  // Return
945  return;
946 }
947 
948 
949 /***********************************************************************//**
950  * @brief Returns index of additional vector column
951  *
952  * @param[in] colname Vector column name.
953  * @return Vector column index.
954  *
955  * Returns index of additional vector column. Returns -1 if the vector column
956  * has not found.
957  ***************************************************************************/
958 int GArf::column_index(const std::string& colname) const
959 {
960  // Initialise index
961  int index(-1);
962 
963  // Search vector column name
964  for (int i = 0; i < columns(); ++i) {
965  if (m_colnames[i] == colname) {
966  index = i;
967  break;
968  }
969  }
970 
971  // Return index
972  return index;
973 }
GArf & operator-=(const GArf &arf)
Subtract Auxiliary Response File.
Definition: GArf.cpp:230
void unit(const std::string &unit)
Set column unit.
const int & ncols(void) const
Return number of columns in table.
Definition: GFitsTable.hpp:134
double & operator[](const int &index)
Return content of spectral bin.
Definition: GArf.hpp:151
#define G_OPERATOR_MINUS
Definition: GArf.cpp:43
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
bool contains(const int &extno) const
Check if HDU exists in FITS file.
Definition: GFits.hpp:282
GArf(void)
Void constructor.
Definition: GArf.cpp:66
void copy_members(const GArf &pha)
Copy class members.
Definition: GArf.cpp:890
GEbounds m_ebounds
Energy boundaries.
Definition: GArf.hpp:122
double & at(const int &index)
Return content of spectral bin.
Definition: GArf.cpp:428
void read(const GFits &fits)
Read Auxiliary Response File.
Definition: GArf.cpp:617
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:163
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition: GEbounds.cpp:301
GFilename m_filename
Filename of origin.
Definition: GArf.hpp:121
bool contains(const std::string &colname) const
Checks the presence of a column in table.
Definition: GFitsTable.cpp:759
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:264
void nrows(const int &nrows)
Set number of rows in column.
#define G_OPERATOR
Definition: GArf.cpp:44
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
void reserve(const int &num)
Reserves space for nodes in node array.
Definition: GNodeArray.hpp:220
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition: GFitsTable.hpp:147
std::vector< std::string > m_colnames
Additional column names.
Definition: GArf.hpp:125
#define G_OPERATOR_PLUS
Definition: GArf.cpp:42
Gammalib tools definition.
FITS table float column class interface definition.
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 save(const GFilename &filename, const bool &clobber=false) const
Save Auxiliary Response File.
Definition: GArf.cpp:589
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:587
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:80
GFitsHeader m_header
FITS header cards.
Definition: GArf.hpp:127
void append(const std::string &name, const std::vector< double > &column)
Append additional column to spectrum.
Definition: GArf.cpp:522
Auxiliary Response File class.
Definition: GArf.hpp:54
#define G_APPEND
Definition: GArf.cpp:48
GArf & operator/=(const double &scale)
Divide Auxiliary Response File values.
Definition: GArf.cpp:277
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1293
GArf & operator=(const GArf &arf)
Assignment operator.
Definition: GArf.cpp:161
int column_index(const std::string &colname) const
Returns index of additional vector column.
Definition: GArf.cpp:958
const std::string extname_arf
Definition: GArf.hpp:45
void free_members(void)
Delete class members.
Definition: GArf.cpp:909
GArf & operator+=(const GArf &arf)
Add Auxiliary Response File.
Definition: GArf.cpp:196
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:279
std::vector< std::vector< double > > m_coldata
Additional column data.
Definition: GArf.hpp:126
Energy boundaries container class.
Definition: GEbounds.hpp:60
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:862
void init_members(void)
Initialise class members.
Definition: GArf.cpp:869
GArf & operator*=(const double &scale)
Scale Auxiliary Response File values.
Definition: GArf.cpp:257
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
Filename class.
Definition: GFilename.hpp:62
std::string print(const GChatter &chatter=NORMAL) const
Print Auxiliary Response File.
Definition: GArf.cpp:834
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:187
Abstract interface for FITS table column.
XSPEC Auxiliary Response File class definition.
virtual ~GArf(void)
Destructor.
Definition: GArf.cpp:139
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
void set_logetrue(void)
Set true energy node array.
Definition: GArf.cpp:919
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
int size(void) const
Return number of cards in header.
GChatter
Definition: GTypemaps.hpp:33
GEnergy elogmean(const int &index) const
Returns logarithmic mean energy for a given energy interval.
Definition: GEbounds.cpp:1126
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:235
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:268
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:249
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
void write(GFits &fits) const
Write Auxiliary Response File.
Definition: GArf.cpp:740
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
#define G_OPERATOR2
Definition: GArf.cpp:45
virtual double real(const int &row, const int &inx=0) const =0
double keV(void) const
Return energy in keV.
Definition: GEnergy.cpp:306
int columns(void) const
Return number of additional columns.
Definition: GArf.hpp:223
FITS binary table class.
void clear(void)
Clear object.
Definition: GArf.cpp:395
#define G_AT2
Definition: GArf.cpp:47
Exception handler interface definition.
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition: GFits.cpp:678
FITS binary table class definition.
std::string tolower(const std::string &s)
Convert string to lower case.
Definition: GTools.cpp:955
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
FITS table short integer column class interface definition.
FITS table float column.
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:199
double & operator()(const int &index, const int &col)
Return content of additional columns.
Definition: GArf.hpp:180
void clear(void)
Clear header.
#define G_AT1
Definition: GArf.cpp:46
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
GNodeArray m_logetrue
Log10 energies in TeV.
Definition: GArf.hpp:123
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
int size(void) const
Return number of spectral bins.
Definition: GArf.hpp:209
void load(const GFilename &filename)
Load Auxiliary Response File.
Definition: GArf.cpp:552
GArf * clone(void) const
Clone object.
Definition: GArf.cpp:411
Filename class interface definition.
const GEbounds & ebounds(void) const
Return energy boundaries.
Definition: GArf.hpp:237
std::vector< double > m_specresp
Spectral response.
Definition: GArf.hpp:124
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:489
FITS table abstract base class interface definition.