GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTACubeBackground.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTACubeBackground.cpp - CTA cube background class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2015-2020 by Michael Mayer *
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 GCTACubeBackground.cpp
23  * @brief CTA cube background class implementation
24  * @author Michael Mayer
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GException.hpp"
33 #include "GMath.hpp"
34 #include "GRan.hpp"
35 #include "GLog.hpp"
36 #include "GFits.hpp"
37 #include "GFitsBinTable.hpp"
38 #include "GObservations.hpp"
39 #include "GSkyRegionCircle.hpp"
40 #include "GCTAEventCube.hpp"
41 #include "GCTAObservation.hpp"
42 #include "GCTARoi.hpp"
43 #include "GCTAInstDir.hpp"
44 #include "GCTACubeBackground.hpp"
45 
46 /* __ Method name definitions ____________________________________________ */
47 #define G_READ "GCTACubeBackground::read(GFits&)"
48 #define G_MC "GCTACubeBackground::mc(GEnergy&, GTime&, GRan&)"
49 #define G_FILL "GCTACubeBackground::fill(GObservations&, GLog*)"
50 
51 /* __ Macros _____________________________________________________________ */
52 
53 /* __ Coding definitions _________________________________________________ */
54 #define G_LOG_INTERPOLATION //!< Energy interpolate log10(background rate)
55 
56 /* __ Debug definitions __________________________________________________ */
57 //#define G_DEBUG_MC_INIT
58 //#define G_DEBUG_CACHE
59 
60 /* __ Constants __________________________________________________________ */
61 
62 
63 /*==========================================================================
64  = =
65  = Constructors/destructors =
66  = =
67  ==========================================================================*/
68 
69 /***********************************************************************//**
70  * @brief Void constructor
71  ***************************************************************************/
73 {
74  // Initialise class members
75  init_members();
76 
77  // Return
78  return;
79 }
80 
81 
82 /***********************************************************************//**
83  * @brief File constructor
84  *
85  * @param[in] filename FITS file name.
86  *
87  * Construct instance by loading the background information from a FITS file.
88  ***************************************************************************/
90 {
91  // Initialise class members
92  init_members();
93 
94  // Load background from file
95  load(filename);
96 
97  // Return
98  return;
99 }
100 
101 
102 /***********************************************************************//**
103  * @brief Event cube constructor
104  *
105  * @param[in] cube Event cube.
106  *
107  * Construct background cube using the same binning and sky projection that
108  * is used for the event cube.
109  ***************************************************************************/
111 {
112  // Initialise class members
113  init_members();
114 
115  // Store energy boundaries
116  m_ebounds = cube.ebounds();
117 
118  // Set GNodeArray used for interpolation
119  set_eng_axis();
120 
121  // Set background cube to event cube
122  m_cube = cube.counts();
124 
125  // Set all background cube pixels to zero as we want to have a clean map
126  // upon construction
127  m_cube = 0.0;
128 
129  // Return
130  return;
131 
132 }
133 
134 
135 /***********************************************************************//**
136  * @brief Copy constructor
137  *
138  * @param[in] bgd Background.
139  ***************************************************************************/
141 {
142  // Initialise class members
143  init_members();
144 
145  // Copy members
146  copy_members(bgd);
147 
148  // Return
149  return;
150 }
151 
152 
153 /***********************************************************************//**
154  * @brief Background cube constructor
155  *
156  * @param[in] wcs World Coordinate System.
157  * @param[in] coords Coordinate System (CEL or GAL).
158  * @param[in] x X coordinate of sky map centre (deg).
159  * @param[in] y Y coordinate of sky map centre (deg).
160  * @param[in] dx Pixel size in x direction at centre (deg/pixel).
161  * @param[in] dy Pixel size in y direction at centre (deg/pixel).
162  * @param[in] nx Number of pixels in x direction.
163  * @param[in] ny Number of pixels in y direction.
164  * @param[in] ebounds Energy boundaries.
165  *
166  * Constructs a background cube by specifying the sky map grid and the
167  * energies.
168  ***************************************************************************/
170  const std::string& coords,
171  const double& x,
172  const double& y,
173  const double& dx,
174  const double& dy,
175  const int& nx,
176  const int& ny,
177  const GEbounds& ebounds)
178 {
179  // Initialise class members
180  init_members();
181 
182  // Store energy boundaries
183  m_ebounds = ebounds;
184 
185  // Set GNodeArray used for interpolation
186  set_eng_axis();
187 
188  // Create sky map
189  m_cube = GSkyMap(wcs, coords, x, y, dx, dy, nx, ny, m_ebounds.size());
190 
191  // Return
192  return;
193 }
194 
195 
196 /***********************************************************************//**
197  * @brief Destructor
198  ***************************************************************************/
200 {
201  // Free members
202  free_members();
203 
204  // Return
205  return;
206 }
207 
208 
209 /*==========================================================================
210  = =
211  = Operators =
212  = =
213  ==========================================================================*/
214 
215 /***********************************************************************//**
216  * @brief Assignment operator
217  *
218  * @param[in] bgd Background cube.
219  * @return Background cube.
220  ***************************************************************************/
222 {
223  // Execute only if object is not identical
224  if (this != &bgd) {
225 
226  // Free members
227  free_members();
228 
229  // Initialise private members
230  init_members();
231 
232  // Copy members
233  copy_members(bgd);
234 
235  } // endif: object was not identical
236 
237  // Return this object
238  return *this;
239 }
240 
241 
242 /***********************************************************************//**
243  * @brief Return background rate in units of events MeV\f$^{-1}\f$
244  * s\f$^{-1}\f$ sr\f$^{-1}\f$
245  *
246  * @param[in] dir Reconstructed event position where rate should be computed
247  * @param[in] energy Energy at which the rate should be computed
248  * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
249  *
250  * Returns the background rate for a giben instrument direction and energy.
251  * The rate is given per ontime. The method assures that the background rate
252  * never becomes negative.
253  *
254  * The method interpolates logarithmically in the energy direction.
255  ***************************************************************************/
257  const GEnergy& energy) const
258 {
259  // Initialise background rate
260  double background = 0.0;
261 
262  // Set indices and weighting factors for interpolation
263  update(energy.log10TeV());
264 
265  // If left weight is close to 1, use left node
266  if (std::abs(m_wgt_left-1.0) < 1.0e-6) {
267  background = m_cube(dir.dir(), m_inx_left);
268  }
269 
270  // ... otherwise, if right weight is close to 1, use right node
271  else if (std::abs(m_wgt_right-1.0) < 1.0e-6) {
272  background = m_cube(dir.dir(), m_inx_right);
273  }
274 
275  // ... otherwise perform interpolation
276  else {
277  double background_left = m_cube(dir.dir(), m_inx_left);
278  double background_right = m_cube(dir.dir(), m_inx_right);
279  if (background_left > 0.0 && background_right > 0.0) {
280  background = std::exp(m_wgt_left * std::log(background_left) +
281  m_wgt_right * std::log(background_right));
282  }
283  }
284 
285  // Make sure that background rate does not become negative
286  if (background < 0.0) {
287  background = 0.0;
288  }
289 
290  // Return background rate
291  return background;
292 }
293 
294 
295 /*==========================================================================
296  = =
297  = Public methods =
298  = =
299  ==========================================================================*/
300 
301 /***********************************************************************//**
302  * @brief Clear background
303  *
304  * This method properly resets the background to an initial state.
305  ***************************************************************************/
307 {
308  // Free class members (base and derived classes, derived class first)
309  free_members();
310 
311  // Initialise members
312  init_members();
313 
314  // Return
315  return;
316 }
317 
318 
319 /***********************************************************************//**
320  * @brief Clone background
321  *
322  * @return Pointer to deep copy of background.
323  ***************************************************************************/
325 {
326  return new GCTACubeBackground(*this);
327 }
328 
329 
330 /***********************************************************************//**
331  * @brief Fill background cube from observation container
332  *
333  * @param[in] obs Observation container.
334  * @param[in] log Pointer to logger (optional).
335  *
336  * @exception GException::invalid_value
337  * No event list found in CTA observations.
338  *
339  * Set the background cube by computing the ontime weighted background rate
340  * for all CTA observations in an observation container. The cube pixel
341  * values are computed as the sum over the background rates.
342  ***************************************************************************/
344 {
345  // Set energy margin
346  const GEnergy energy_margin(0.01, "GeV");
347 
348  // Clear background cube
349  m_cube = 0.0;
350 
351  // Set dummy GTI needed to generate an event cube. It is not important
352  // what the actually value is since it will be overwritten later in any
353  // case, but it's important that there is one time slice
354  GGti gti(GTime(0.0), GTime(1.0));
355 
356  // Initialise event cubes for model evaluate and storage
357  GCTAEventCube eventcube = GCTAEventCube(m_cube, m_ebounds, gti);
358 
359  // Initialise total ontime
360  double total_ontime = 0.0;
361 
362  // Loop over all observations in container
363  for (int i = 0; i < obs.size(); ++i) {
364 
365  // Get observation and continue only if it is a CTA observation
366  const GCTAObservation* cta =
367  dynamic_cast<const GCTAObservation*>(obs[i]);
368 
369  // Skip observation if it's not CTA
370  if (cta == NULL) {
371  if (log != NULL) {
372  *log << "Skipping ";
373  *log << obs[i]->instrument();
374  *log << " observation ";
375  *log << "\"" << obs[i]->name() << "\"";
376  *log << " (id=" << obs[i]->id() << ")" << std::endl;
377  }
378  continue;
379  }
380 
381  // Skip observation if we have a binned observation
382  if (cta->eventtype() == "CountsCube") {
383  if (log != NULL) {
384  *log << "Skipping binned ";
385  *log << cta->instrument();
386  *log << " observation ";
387  *log << "\"" << cta->name() << "\"";
388  *log << " (id=" << cta->id() << ")" << std::endl;
389  }
390  continue;
391  }
392 
393  // Get observation ontime
394  double ontime = cta->ontime();
395 
396  // Skip observation if ontime is zero
397  if (ontime == 0.0) {
398  if (log != NULL) {
399  *log << "Skipping unbinned ";
400  *log << cta->instrument();
401  *log << " observation ";
402  *log << "\"" << cta->name() << "\"";
403  *log << " (id=" << cta->id() << ") due to zero ontime";
404  *log << std::endl;
405  }
406  continue;
407  }
408 
409  // Extract region of interest from CTA observation
410  GCTARoi roi = cta->roi();
411 
412  // Extract energy boundaries from CTA observation
413  GEbounds obs_ebounds = cta->ebounds();
414 
415  // Check for RoI sanity
416  if (!roi.is_valid()) {
417  std::string msg = "No RoI information found in input observation "
418  "\""+cta->name()+"\". Run ctselect to specify "
419  "an RoI for this observation";
420  throw GException::invalid_value(G_FILL, msg);
421  }
422 
423  // Make sure the observation overlaps with the cube
424  GSkyRegionCircle obs_reg(roi.centre().dir(), roi.radius());
425  if (!m_cube.overlaps(obs_reg)) {
426  if (log != NULL) {
427  *log << "Skipping unbinned ";
428  *log << cta->instrument();
429  *log << " observation ";
430  *log << "\"" << cta->name() << "\"";
431  *log << " (id=" << cta->id() << ") since it does not overlap ";
432  *log << "with the background cube.";
433  *log << std::endl;
434  }
435  continue;
436  } // endif: m_cube overlaps observation
437 
438  // Announce observation usage
439  if (log != NULL) {
440  *log << "Including ";
441  *log << cta->instrument();
442  *log << " observation \"" << cta->name();
443  *log << "\" (id=" << cta->id() << ")";
444  *log << " in background cube computation." << std::endl;
445  }
446 
447  // Set GTI of actual observations as the GTI of the event cube
448  eventcube.gti(cta->gti());
449 
450  // Compute number of pixels and energy layers
451  int npix = eventcube.npix();
452  int nebins = m_ebounds.size();
453 
454  // Loop over all spatial pixels
455  for (int ipix = 0; ipix < npix; ++ipix) {
456 
457  // Get event bin in first layer
458  GCTAEventBin* bin = eventcube[ipix];
459 
460  // Skip if bin is not contained in RoI
461  if (!roi.contains(*bin)) {
462  continue;
463  }
464 
465  // Get CTA direction for this bin, including sky direction and
466  // DETX/DETY coordinates
467  GCTAInstDir dir = cta->pointing().instdir(bin->dir().dir());
468 
469  // Loop over all energy bins of the cube
470  for (int iebin = 0, ibin = ipix; iebin < nebins; ++iebin, ibin += npix) {
471 
472  // Compute background model value if the energy bin is fully
473  // contained in the observation interval. Add a small margin
474  // here to the energy bin boundaries to take provision for
475  // rounding errors.
476  if (obs_ebounds.contains(m_ebounds.emin(iebin) + energy_margin,
477  m_ebounds.emax(iebin) - energy_margin)) {
478 
479  // Set energy bin
480  GCTAEventBin* bin = eventcube[ibin];
481  bin->dir(dir); // Set instrument direction
482  bin->energy(m_ebounds.elogmean(iebin));
483 
484  // Compute model
485  double model = obs.models().eval(*bin, *cta) * ontime;
486 
487  // Add existing number of counts
488  model += bin->counts();
489 
490  // Store cumulated value (units: counts/(MeV sr))
491  bin->counts(model);
492 
493  } // endif: added background model value to cube
494 
495  } // endfor: looped over energy bins
496 
497  } // endfor: looped over spatial pixels
498 
499  // Accumulate ontime
500  total_ontime += ontime;
501 
502  } // endfor: looped over all observations
503 
504  // Re-normalize cube to get units of counts/(MeV s sr)
505  if (total_ontime > 0.0) {
506 
507  // Loop over all bins in background cube and divide the content
508  // by the total ontime.
509  for (int i = 0; i < eventcube.size(); ++i) {
510  GCTAEventBin* bin = eventcube[i];
511  double rate = bin->counts() / total_ontime;
512  bin->counts(rate);
513  }
514 
515  // Set background cube values from event cube
516  m_cube = eventcube.counts();
517 
518  } // endif: ontime was positive
519 
520  // Return
521  return;
522 
523 }
524 
525 
526 /***********************************************************************//**
527  * @brief Return spatially integrated background rate in units of
528  * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
529  *
530  * @param[in] logE Logarithm (base 10) of energy in TeV
531  * @return Spatially integrated background rate
532  * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
533  *
534  * Spatially integrates the background cube at a given energy. This method
535  * performs an interpolation between the energy maps. The integration is
536  * performed by summing over all bin contents multiplied by the solid angle
537  * of the bin.
538  ***************************************************************************/
539 double GCTACubeBackground::integral(const double& logE) const
540 {
541  // Update interpolation cache
542  update(logE);
543 
544  // Initialise result
545  double result = 0.0;
546 
547  // Loop over all map pixels
548  for (int i = 0; i < m_cube.npix(); ++i) {
549 
550  // Get bin value
551  double value = m_wgt_left * m_cube(i, m_inx_left) +
553 
554  // Sum bin contents
555  result += value * m_cube.solidangle(i);
556 
557  }
558 
559  // Return result
560  return result;
561 }
562 
563 
564 /***********************************************************************//**
565  * @brief Read background cube from FITS object
566  *
567  * @param[in] fits FITS object.
568  *
569  * Read the background cube from a FITS object.
570  ***************************************************************************/
572 {
573  // Clear object
574  clear();
575 
576  // Get HDUs
577  const GFitsImage& hdu_bgdcube = *fits.image("Primary");
578  const GFitsTable& hdu_ebounds = *fits.table(gammalib::extname_ebounds);
579 
580  // Read cube
581  m_cube.read(hdu_bgdcube);
582 
583  // Read energy boundaries
584  m_ebounds.read(hdu_ebounds);
585 
586  // Set energy node array
587  set_eng_axis();
588 
589  // Return
590  return;
591 }
592 
593 
594 /***********************************************************************//**
595  * @brief Write CTA background cube into FITS object.
596  *
597  * @param[in] fits FITS file.
598  ***************************************************************************/
600 {
601  // Write cube
602  m_cube.write(fits);
603 
604  // Get last HDU and write attributes
605  if (fits.size() > 0) {
606  GFitsHDU& hdu = *fits[fits.size()-1];
607  hdu.card("BUNIT", "MeV**(-1) s**(-1) sr**(-1)",
608  "Unit of background cube");
609  }
610 
611  // Write energy boundaries
612  m_ebounds.write(fits);
613 
614  // Return
615  return;
616 }
617 
618 
619 /***********************************************************************//**
620  * @brief Load background cube from FITS file
621  *
622  * @param[in] filename FITS file.
623  *
624  * Loads the background cube from a FITS file.
625  ***************************************************************************/
626 void GCTACubeBackground::load(const GFilename& filename)
627 {
628  // Put into OpenMP criticial zone
629  #pragma omp critical(GCTACubeBackground_load)
630  {
631 
632  // Open FITS file
633  GFits fits(filename);
634 
635  // Read background cube
636  read(fits);
637 
638  // Close FITS file
639  fits.close();
640 
641  } // end of OpenMP critical zone
642 
643  // Store filename
645 
646  // Return
647  return;
648 
649 }
650 
651 
652 /***********************************************************************//**
653  * @brief Save background cube into FITS file
654  *
655  * @param[in] filename background cube FITS file name.
656  * @param[in] clobber Overwrite existing file?
657  *
658  * Save the background cube into a FITS file.
659  ***************************************************************************/
660 void GCTACubeBackground::save(const GFilename& filename,
661  const bool& clobber) const
662 {
663  // Put into OpenMP criticial zone
664  #pragma omp critical(GCTACubeBackground_save)
665  {
666 
667  // Open or create FITS file
668  GFits fits;
669 
670  // Write background cube
671  write(fits);
672 
673  // Save FITS file
674  fits.saveto(filename, clobber);
675 
676  // Close Edisp file
677  fits.close();
678 
679  } // end of OpenMP critical zone
680 
681  // Store filename
683 
684  // Return
685  return;
686 }
687 
688 
689 /***********************************************************************//**
690  * @brief Print background information
691  *
692  * @param[in] chatter Chattiness.
693  * @return String containing background information.
694  ***************************************************************************/
695 std::string GCTACubeBackground::print(const GChatter& chatter) const
696 {
697  // Initialise result string
698  std::string result;
699 
700  // Continue only if chatter is not silent
701  if (chatter != SILENT) {
702 
703  // Append header
704  result.append("=== GCTACubeBackground ===");
705 
706  // Append parameters
707  result.append("\n"+gammalib::parformat("Filename")+m_filename);
708 
709  // Append energies
710  if (m_ebounds.size() > 0) {
711  result.append("\n"+m_ebounds.print(chatter));
712  }
713  else {
714  result.append("\n"+gammalib::parformat("Energy boundaries") +
715  "Not defined");
716  }
717 
718  // Append skymap definition
719  result.append("\n"+m_cube.print(chatter));
720 
721  } // endif: chatter was not silent
722 
723  // Return result
724  return result;
725 }
726 
727 
728 /*==========================================================================
729  = =
730  = Private methods =
731  = =
732  ==========================================================================*/
733 
734 /***********************************************************************//**
735  * @brief Initialise class members
736  ***************************************************************************/
738 {
739  // Initialise members
740  m_filename.clear();
741  m_cube.clear();
742  m_ebounds.clear();
743  m_elogmeans.clear();
744 
745  // Initialise cache
746  m_inx_left = 0;
747  m_inx_right = 0;
748  m_wgt_left = 0.0;
749  m_wgt_right = 0.0;
750 
751  // Return
752  return;
753 }
754 
755 
756 /***********************************************************************//**
757  * @brief Copy class members
758  *
759  * @param[in] bgd Background.
760  ***************************************************************************/
762 {
763  // Copy members
764  m_filename = bgd.m_filename;
765  m_cube = bgd.m_cube;
766  m_ebounds = bgd.m_ebounds;
767  m_elogmeans = bgd.m_elogmeans;
768 
769  // Copy cache
770  m_inx_left = bgd.m_inx_left;
771  m_inx_right = bgd.m_inx_right;
772  m_wgt_left = bgd.m_wgt_left;
773  m_wgt_right = bgd.m_wgt_right;
774 
775  // Return
776  return;
777 }
778 
779 
780 /***********************************************************************//**
781  * @brief Delete class members
782  ***************************************************************************/
784 {
785  // Return
786  return;
787 }
788 
789 
790 /***********************************************************************//**
791  * @brief Set nodes for a logarithmic (base 10) energy axis
792  *
793  *
794  * Set axis nodes so that each node is the logarithm of the energy values.
795  ***************************************************************************/
797 {
798  // Get number of bins
799  int bins = m_ebounds.size();
800 
801  // Clear node array
802  m_elogmeans.clear();
803 
804  // Compute nodes
805  for (int i = 0; i < bins; ++i) {
806 
807  // Get logE/TeV
809 
810  } // endfor: looped over energy bins
811 
812  // Return
813  return;
814 }
815 
816 
817 /***********************************************************************//**
818  * @brief Update 1D cache
819  *
820  * @param[in] logE Log10 energy in TeV.
821  *
822  * Updates the 1D interpolation cache. The interpolation cache is composed
823  * of two indices and weights that define 2 data values of the 2D skymap
824  * that are used for linear interpolation.
825  *
826  * @todo Write down formula
827  ***************************************************************************/
828 void GCTACubeBackground::update(const double& logE) const
829 {
830  // Set value for node array
831  m_elogmeans.set_value(logE);
832 
833  // Set indices and weighting factors for interpolation
838 
839  // Return
840  return;
841 }
Sky map class.
Definition: GSkyMap.hpp:89
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
void init_members(void)
Initialise class members.
GEbounds m_ebounds
Energy boundaries of background cube.
void free_members(void)
Delete class members.
Random number generator class definition.
virtual const GCTAInstDir & dir(void) const
Return instrument direction of event bin.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
const GFilename & filename(void) const
Return background cube filename.
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
bool overlaps(const GSkyRegion &region) const
Checks whether a region overlaps with this map.
Definition: GSkyMap.cpp:2102
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:163
GCTAEventBin class interface definition.
CTA cube background class.
GCTACubeBackground(void)
Void constructor.
GEbounds ebounds(void) const
Get energy boundaries.
bool contains(const GEnergy &eng) const
Checks whether energy boundaries contain energy.
Definition: GEbounds.cpp:1180
const double & radius(void) const
Returns radius of region of interest in degrees.
Definition: GCTARoi.hpp:125
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:264
void gti(const GGti &gti)
Set Good Time Intervals.
Definition: GEvents.cpp:154
void models(const GModels &models)
Set model container.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
void counts(const GSkyMap &counts)
Set event cube counts from sky map.
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
std::string print(const GChatter &chatter=NORMAL) const
Print background information.
Interface for the CTA region of interest class.
Definition: GCTARoi.hpp:49
void read(const GFitsTable &table)
Read energy boundaries from FITS table.
Definition: GEbounds.cpp:761
Time class.
Definition: GTime.hpp:55
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
virtual const GEnergy & energy(void) const
Return energy of event bin.
FITS file class interface definition.
Interface for the circular sky region class.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:587
void pointing(const GCTAPointing &pointing)
Set CTA pointing.
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2664
void dir(const GSkyDir &dir)
Set sky direction.
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:816
void update(const double &logE) const
Update 1D cache.
void id(const std::string &id)
Set observation identifier.
void write(GFits &file) const
Write CTA background cube into FITS object.
int m_inx_right
Index of right node.
virtual ~GCTACubeBackground(void)
Destructor.
GCTACubeBackground * clone(void) const
Clone background.
GFilename m_filename
Name of background response file.
CTA cube background class definition.
virtual double ontime(void) const
Return ontime.
double m_wgt_left
Weight of left node.
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1293
std::string print(const GChatter &chatter=NORMAL) const
Print sky map.
Definition: GSkyMap.cpp:2787
void read(const GFits &fits)
Read background cube from FITS object.
Information logger interface definition.
Definition: GLog.hpp:62
void set_eng_axis(void)
Set nodes for a logarithmic (base 10) energy axis.
CTA event bin container class.
GNodeArray m_elogmeans
Mean energy for the background cube.
void clear(void)
Clear background.
const GEbounds & ebounds(void) const
Return energy boundaries.
int npix(void) const
Return number of pixels in one energy bins of the event cube.
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:279
Energy boundaries container class.
Definition: GEbounds.hpp:60
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:389
GCTARoi roi(void) const
Get Region of Interest.
GFitsHDU * write(GFits &file, const std::string &extname="") const
Write sky map into FITS file.
Definition: GSkyMap.cpp:2697
Filename class.
Definition: GFilename.hpp:62
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:187
GGti gti(void) const
Get Good Time Intervals.
const GCTAInstDir & centre(void) const
Returns region of interest centre.
Definition: GCTARoi.hpp:111
virtual std::string instrument(void) const
Return instrument name.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
const std::string extname_ebounds
Definition: GEbounds.hpp:44
CTA instrument direction class interface definition.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
virtual bool contains(const GEvent &event) const
Check if region of interest contains an event.
Definition: GCTARoi.cpp:215
std::string print(const GChatter &chatter=NORMAL) const
Print energy boundaries.
Definition: GEbounds.cpp:1232
GChatter
Definition: GTypemaps.hpp:33
int m_inx_left
Index of left node.
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
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1100
Good Time Interval class.
Definition: GGti.hpp:62
double integral(const double &logE) const
Return spatially integrated background rate in units of events MeV s .
int size(void) const
Return number of observations in container.
void copy_members(const GCTACubeBackground &bgd)
Copy class members.
GCTACubeBackground & operator=(const GCTACubeBackground &bgd)
Assignment operator.
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:268
CTA region of interest class interface definition.
CTA observation class interface definition.
Observation container class.
#define G_FILL
void load(const GFilename &filename)
Load background cube from FITS file.
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:249
void name(const std::string &name)
Set observation name.
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:368
Information logger class definition.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
void fill(const GObservations &obs, GLog *log=NULL)
Fill background cube from observation container.
int size(void) const
Return number of HDUs in FITS file.
Definition: GFits.hpp:237
void save(const GFilename &filename, const bool &clobber=false) const
Save background cube into FITS file.
CTA event bin container class interface definition.
std::string eventtype(void) const
Return event type string.
bool is_valid(void) const
Checks if RoI is valid.
Definition: GCTARoi.hpp:152
Exception handler interface definition.
FITS binary table class definition.
double m_wgt_right
Weight of right node.
double operator()(const GCTAInstDir &dir, const GEnergy &energy) const
Return background rate in units of events MeV s sr .
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
CTA instrument direction class.
Definition: GCTAInstDir.hpp:63
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:199
CTA observation class.
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
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:345
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
Observations container class interface definition.
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
virtual double counts(void) const
Return number of counts in event bin.
Circular sky region class interface definition.
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
GSkyMap m_cube
Background cube.
virtual int size(void) const
Return number of bins in event cube.