GammaLib  1.7.0.dev
 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-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 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 Bin index [0,...,size()-1].
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, index, 0, size()-1);
433  }
434 
435  // Return reference
436  return (m_specresp[index]);
437 }
438 
439 
440 /***********************************************************************//**
441  * @brief Return content of spectral bin (const version)
442  *
443  * @param[in] index Bin index [0,...,size()-1].
444  *
445  * @exception GException::out_of_range
446  * Spectral bin index is out of range.
447  *
448  * Returns reference to content of spectral bin with specified @p index.
449  ***************************************************************************/
450 const double& GArf::at(const int& index) const
451 {
452  // Throw an exception if index is out of range
453  if (index < 0 || index >= size()) {
454  throw GException::out_of_range(G_AT1, index, 0, size()-1);
455  }
456 
457  // Return reference
458  return (m_specresp[index]);
459 }
460 
461 
462 /***********************************************************************//**
463  * @brief Return content of additional columns
464  *
465  * @param[in] index Bin index [0,...,size()-1].
466  * @param[in] col Columns index [0,...,columns()-1].
467  *
468  * @exception GException::out_of_range
469  * Spectral bin or column index is out of range.
470  *
471  * Returns reference to content of additional columns.
472  ***************************************************************************/
473 double& GArf::at(const int& index, const int& col)
474 {
475  // Throw an exception if bin or column index is out of range
476  if (index < 0 || index >= size()) {
477  throw GException::out_of_range(G_AT2, "Bin index", index, size());
478  }
479  if (col < 0 || col >= columns()) {
480  throw GException::out_of_range(G_AT2, "Column index", col, columns());
481  }
482 
483  // Return reference
484  return (m_coldata[col][index]);
485 }
486 
487 
488 /***********************************************************************//**
489  * @brief Return content of additional columns (const version)
490  *
491  * @param[in] index Bin index [0,...,size()-1].
492  * @param[in] col Columns index [0,...,columns()-1].
493  *
494  * @exception GException::out_of_range
495  * Spectral bin or column index is out of range.
496  *
497  * Returns reference to content of additional columns.
498  ***************************************************************************/
499 const double& GArf::at(const int& index, const int& col) const
500 {
501  // Throw an exception if bin or column index is out of range
502  if (index < 0 || index >= size()) {
503  throw GException::out_of_range(G_AT2, "Bin index", index, size());
504  }
505  if (col < 0 || col >= columns()) {
506  throw GException::out_of_range(G_AT2, "Column index", col, columns());
507  }
508 
509  // Return reference
510  return (m_coldata[col][index]);
511 }
512 
513 
514 /***********************************************************************//**
515  * @brief Append additional column to spectrum
516  *
517  * @param[in] name Additional column name.
518  * @param[in] column Additional column data.
519  ***************************************************************************/
520 void GArf::append(const std::string& name, const std::vector<double>& column)
521 {
522  // Throw an exception if the number of elements in the column does not
523  // correspond to the size of the spectrum
524  if (column.size() != size()) {
525  std::string msg = "Size of column "+gammalib::str(column.size())+
526  " is incompatible with size of spectrum "+
527  gammalib::str(size())+".";
529  }
530 
531  // Append column name
532  m_colnames.push_back(name);
533 
534  // Append column data
535  m_coldata.push_back(column);
536 
537  // Return
538  return;
539 }
540 
541 
542 /***********************************************************************//**
543  * @brief Load Auxiliary Response File
544  *
545  * @param[in] filename File name.
546  *
547  * Loads the Auxiliary Response File from the `SPECRESP` extension of the
548  * FITS file.
549  ***************************************************************************/
550 void GArf::load(const GFilename& filename)
551 {
552  // Clear response
553  clear();
554 
555  // Open FITS file (without extension name as the user is not allowed
556  // to modify the extension names)
557  GFits fits(filename.url());
558 
559  // Read ARF data
560  read(fits);
561 
562  // Close FITS file
563  fits.close();
564 
565  // Store filename
566  m_filename = filename.url();
567 
568  // Return
569  return;
570 }
571 
572 
573 /***********************************************************************//**
574  * @brief Save Auxiliary Response File
575  *
576  * @param[in] filename File name.
577  * @param[in] clobber Overwrite existing file?
578  *
579  * Saves the Auxiliary Response File into a FITS file. If a file with the
580  * given @p filename does not yet exist it will be created. If the file
581  * exists it can be overwritten if the @p clobber flag is set to `true`.
582  * Otherwise an exception is thrown.
583  *
584  * The method will save the `SPECRESP` binary FITS table into the FITS file
585  * that contains the values of the Auxiliary Response File.
586  ***************************************************************************/
587 void GArf::save(const GFilename& filename, const bool& clobber) const
588 {
589  // Create FITS file
590  GFits fits;
591 
592  // Write ARF into file
593  write(fits);
594 
595  // Save to file (without extension name since the requested extension
596  // may not yet exist in the file)
597  fits.saveto(filename.url(), clobber);
598 
599  // Store filename
600  m_filename = filename.url();
601 
602  // Return
603  return;
604 }
605 
606 
607 /***********************************************************************//**
608  * @brief Read Auxiliary Response File
609  *
610  * @param[in] fits FITS file.
611  *
612  * Loads the Auxiliary Response File from the `SPECRESP` extension of the
613  * FITS file.
614  ***************************************************************************/
615 void GArf::read(const GFits& fits)
616 {
617  // Clear response
618  clear();
619 
620  // Get ARF table
621  const GFitsTable& table = *fits.table(gammalib::extname_arf);
622 
623  // Read ARF data
624  read(table);
625 
626  // Return
627  return;
628 }
629 
630 
631 /***********************************************************************//**
632  * @brief Read Auxiliary Response File
633  *
634  * @param[in] table ARF FITS table.
635  *
636  * Reads the Auxiliary Response File from a FITS table. The true energy
637  * boundaries are expected in the `ENERG_LO` and `ENERG_HI` columns, the
638  * response information is expected in the `SPECRESP` column.
639  *
640  * The method will analyze the unit of the `SPECRESP` column, and if either
641  * `m**2`, `m^2` or `m2` are encountered, multiply the values of the column
642  * by \f$10^4\f$ to convert the response into units of \f$cm^2\f$. Units of
643  * the `ENERG_LO` and `ENERG_HI` columns are also interpreted for conversion.
644  *
645  * See
646  * http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html#tth_sEc4
647  * for details about the Auxiliary Response File format.
648  ***************************************************************************/
649 void GArf::read(const GFitsTable& table)
650 {
651  // Clear members
652  clear();
653 
654  // Get pointer to data columns
655  const GFitsTableCol* energy_lo = table["ENERG_LO"];
656  const GFitsTableCol* energy_hi = table["ENERG_HI"];
657  const GFitsTableCol* specresp = table["SPECRESP"];
658 
659  // Determine effective area conversion factor. Internal
660  // units are cm^2
661  std::string u_specresp =
663  double c_specresp = 1.0;
664  if (u_specresp == "m**2" || u_specresp == "m^2" || u_specresp == "m2") {
665  c_specresp = 10000.0;
666  }
667 
668  // Extract number of energy bins
669  int num = energy_lo->nrows();
670 
671  // Set energy bins
672  for (int i = 0; i < num; ++i) {
673 
674  // Append energy bin
675  GEnergy emin(energy_lo->real(i), energy_lo->unit());
676  GEnergy emax(energy_hi->real(i), energy_hi->unit());
677  m_ebounds.append(emin, emax);
678 
679  // Append effective area value
680  double aeff = specresp->real(i) * c_specresp;
681  m_specresp.push_back(aeff);
682 
683  } // endfor: looped over energy bins
684 
685  // Read any additional columns
686  for (int icol = 0; icol < table.ncols(); ++icol) {
687 
688  // Fall through if the column is a standard column
689  std::string colname(table[icol]->name());
690  if ((colname == "ENERG_LO") ||
691  (colname == "ENERG_HI") ||
692  (colname == "SPECRESP")) {
693  continue;
694  }
695 
696  // Get pointer to column
697  const GFitsTableCol* column = table[icol];
698 
699  // Set column vector
700  std::vector<double> coldata;
701  for (int i = 0; i < num; ++i) {
702  coldata.push_back(column->real(i));
703  }
704 
705  // Append column
706  append(colname, coldata);
707 
708  } // endfor: looped over all additional columns
709 
710  // Set log true energy node array
711  set_logetrue();
712 
713  // Store FITS header
714  m_header = table.header();
715 
716  // Return
717  return;
718 }
719 
720 
721 /***********************************************************************//**
722  * @brief Write Auxiliary Response File
723  *
724  * @param[in] fits FITS file.
725  *
726  * Writes the Auxiliary Response File into the `SPECRESP` extension of the
727  * FITS file. An existing extension with the same name will be removed from
728  * the FITS file before writing.
729  *
730  * The method writes the boundaries of the true energy bins into the
731  * `ENERG_LO` and `ENERG_HI` columns, response information will be written
732  * into the `SPECRESP` column.
733  *
734  * See
735  * http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html#tth_sEc4
736  * for details about the Auxiliary Response File format.
737  ***************************************************************************/
738 void GArf::write(GFits& fits) const
739 {
740  // If the FITS object contains already an extension with the same
741  // name then remove now this extension
742  if (fits.contains(gammalib::extname_arf)) {
744  }
745 
746  // Set column length
747  int length = size();
748 
749  // Continue only if there are bins
750  if (length > 0) {
751 
752  // Create new binary table
753  GFitsBinTable hdu;
754 
755  // Allocate floating point vector columns
756  GFitsTableFloatCol energy_lo("ENERG_LO", length);
757  GFitsTableFloatCol energy_hi("ENERG_HI", length);
758  GFitsTableFloatCol specresp("SPECRESP", length);
759 
760  // Fill columns
761  for (int i = 0; i < length; ++i) {
762  energy_lo(i) = (float)m_ebounds.emin(i).keV();
763  energy_hi(i) = (float)m_ebounds.emax(i).keV();
764  specresp(i) = (float)m_specresp[i];
765  }
766 
767  // Set column units
768  energy_lo.unit("keV");
769  energy_hi.unit("keV");
770  specresp.unit("cm**2");
771 
772  // Set table attributes
774 
775  // Append columns to table
776  hdu.append(energy_lo);
777  hdu.append(energy_hi);
778  hdu.append(specresp);
779 
780  // Append any additional columns
781  for (int icol = 0; icol < columns(); ++icol) {
782 
783  // Allocate floating point vector columns
784  GFitsTableFloatCol column(m_colnames[icol], length);
785 
786  // Fill columns
787  for (int i = 0; i < length; ++i) {
788  column(i) = (float)m_coldata[icol][i];
789  }
790 
791  // Append columns to table
792  hdu.append(column);
793 
794  } // endfor: looped over all additional columns
795 
796  // Write mandatory header keywords
797  hdu.card("TELESCOP", "unknown", "Telescope");
798  hdu.card("INSTRUME", "unknown", "Instrument");
799  hdu.card("FILTER", "none", "Filter");
800  hdu.card("HDUCLASS", "OGIP", "Format conforms to OGIP standard");
801  hdu.card("HDUCLAS1", "RESPONSE", "Extension contains response data");
802  hdu.card("HDUCLAS2", "SPECRESP", "Extension contains an ARF");
803  hdu.card("HDUVERS", "1.1.0", "Version of the file format");
804 
805  // Write additional header keywords (do not overwrite any existing
806  // keywords, except of "TELESCOP", "INSTRUME" and "FILTER")
807  for (int i = 0; i < m_header.size(); ++i) {
808  if ((!hdu.contains(m_header[i].keyname())) ||
809  (m_header[i].keyname() == "TELESCOP") ||
810  (m_header[i].keyname() == "INSTRUME") ||
811  (m_header[i].keyname() == "FILTER")) {
812  hdu.card(m_header[i]);
813  }
814  }
815 
816  // Append HDU to FITS file
817  fits.append(hdu);
818 
819  } // endif: there were data to write
820 
821  // Return
822  return;
823 }
824 
825 
826 /***********************************************************************//**
827  * @brief Print Auxiliary Response File
828  *
829  * @param[in] chatter Chattiness.
830  * @return String containing Auxiliary Response File information.
831  ***************************************************************************/
832 std::string GArf::print(const GChatter& chatter) const
833 {
834  // Initialise result string
835  std::string result;
836 
837  // Continue only if chatter is not silent
838  if (chatter != SILENT) {
839 
840  // Append header
841  result.append("=== GArf ===");
842 
843  // Append energy boundary information
844  result.append("\n"+gammalib::parformat("Number of bins"));
845  result.append(gammalib::str(m_ebounds.size()));
846  result.append("\n"+gammalib::parformat("Energy range"));
847  result.append(m_ebounds.emin().print());
848  result.append(" - ");
849  result.append(m_ebounds.emax().print());
850 
851  } // endif: chatter was not silent
852 
853  // Return result
854  return result;
855 }
856 
857 
858 /*==========================================================================
859  = =
860  = Private methods =
861  = =
862  ==========================================================================*/
863 
864 /***********************************************************************//**
865  * @brief Initialise class members
866  ***************************************************************************/
868 {
869  // Initialise members
870  m_filename.clear();
871  m_ebounds.clear();
872  m_logetrue.clear();
873  m_specresp.clear();
874  m_colnames.clear();
875  m_coldata.clear();
876  m_header.clear();
877 
878  // Return
879  return;
880 }
881 
882 
883 /***********************************************************************//**
884  * @brief Copy class members
885  *
886  * @param[in] arf Auxiliary Response File.
887  ***************************************************************************/
888 void GArf::copy_members(const GArf& arf)
889 {
890  // Copy members
891  m_filename = arf.m_filename;
892  m_ebounds = arf.m_ebounds;
893  m_logetrue = arf.m_logetrue;
894  m_specresp = arf.m_specresp;
895  m_colnames = arf.m_colnames;
896  m_coldata = arf.m_coldata;
897  m_header = arf.m_header;
898 
899  // Return
900  return;
901 }
902 
903 
904 /***********************************************************************//**
905  * @brief Delete class members
906  ***************************************************************************/
908 {
909  // Return
910  return;
911 }
912 
913 
914 /***********************************************************************//**
915  * @brief Set true energy node array
916  ***************************************************************************/
918 {
919  // Clear node array
920  m_logetrue.clear();
921 
922  // Continue only if there are true energies in Arf
923  int netrue = m_ebounds.size();
924  if (netrue > 0) {
925 
926  // Reserve space in node array
927  m_logetrue.reserve(netrue);
928 
929  // Append all log mean energies to node array
930  for (int i = 0; i < netrue; ++i) {
931 
932  // Get log mean of true energy in TeV
933  double logE = m_ebounds.elogmean(i).log10TeV();
934 
935  // Append energy to node array
936  m_logetrue.append(logE);
937 
938  } // endfor: appended log mean energies
939 
940  } // endif: there were true energies in Arf
941 
942  // Return
943  return;
944 }
945 
946 
947 /***********************************************************************//**
948  * @brief Returns index of additional vector column
949  *
950  * @param[in] colname Vector column name.
951  * @return Vector column index.
952  *
953  * Returns index of additional vector column. Returns -1 if the vector column
954  * has not found.
955  ***************************************************************************/
956 int GArf::column_index(const std::string& colname) const
957 {
958  // Initialise index
959  int index(-1);
960 
961  // Search vector column name
962  for (int i = 0; i < columns(); ++i) {
963  if (m_colnames[i] == colname) {
964  index = i;
965  break;
966  }
967  }
968 
969  // Return index
970  return index;
971 }
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:472
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:888
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:615
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition: GEbounds.cpp:302
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:716
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
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:216
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:587
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:73
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:520
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:1267
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:956
const std::string extname_arf
Definition: GArf.hpp:45
void free_members(void)
Delete class members.
Definition: GArf.cpp:907
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:275
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:848
void init_members(void)
Initialise class members.
Definition: GArf.cpp:867
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:832
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
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:1184
void set_logetrue(void)
Set true energy node array.
Definition: GArf.cpp:917
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:1140
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:269
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
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:738
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:665
FITS binary table class definition.
std::string tolower(const std::string &s)
Convert string to lower case.
Definition: GTools.cpp:834
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
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:193
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:1022
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:550
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:413
FITS table abstract base class interface definition.