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