GammaLib  1.7.0.dev
 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-2018 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  // Initialise values for integration over energy bin
470  double f1 = 0.0;
471  double f2 = 0.0;
472 
473  // Loop over all energy bins of the cube
474  for (int iebin = 0, ibin = ipix; iebin < nebins; ++iebin, ibin += npix) {
475 
476  // If this is the first bin then compute f1 at the lower bin edge
477  if (iebin == 0) {
478  GCTAEventBin* bin = eventcube[ibin];
479  bin->dir(dir); // Set instrument direction
480  bin->energy(m_ebounds.emin(iebin));
481  f1 = obs.models().eval(*bin, *cta);
482  }
483 
484  // Compute f2 at the upper bin edge
485  GCTAEventBin* bin = eventcube[ibin];
486  bin->dir(dir); // Set instrument direction
487  bin->energy(m_ebounds.emax(iebin));
488  f2 = obs.models().eval(*bin, *cta);
489 
490  // Add background model value to cube if the energy bin is
491  // fully container in the observation interval and if both f1
492  // and f2 are positive. Add a small margin here to the energy
493  // bin boundaries to take provision for rounding errors.
494  if (obs_ebounds.contains(m_ebounds.emin(iebin) + energy_margin,
495  m_ebounds.emax(iebin) - energy_margin) &&
496  (f1 > 0.0 && f2 > 0.0)) {
497 
498  // Compute background model value by integrating over the
499  // bin and deviding the result by the bin width
500  double x1 = m_ebounds.emin(iebin).MeV();
501  double x2 = m_ebounds.emax(iebin).MeV();
502  double model = gammalib::plaw_integral(x1, f1, x2, f2) / (x2 - x1);
503 
504  // Multiply by ontime to get the correct weighting for each
505  // observation. We divide by the total ontime later to get
506  // the background model in units of counts/(MeV s sr).
507  model *= ontime;
508 
509  // Add existing number of counts
510  model += bin->counts();
511 
512  // Store cumulated value (units: counts/(MeV sr))
513  bin->counts(model);
514 
515  } // endif: added background model value to cube
516 
517  // Set f1 to f2 for next energy layer
518  f1 = f2;
519 
520  } // endfor: looped over energy layers
521 
522  } // endfor: looped over spatial pixels
523 
524  // Accumulate ontime
525  total_ontime += ontime;
526 
527  } // endfor: looped over all observations
528 
529  // Re-normalize cube to get units of counts/(MeV s sr)
530  if (total_ontime > 0.0) {
531 
532  // Loop over all bins in background cube and divide the content
533  // by the total ontime.
534  for (int i = 0; i < eventcube.size(); ++i) {
535  GCTAEventBin* bin = eventcube[i];
536  double rate = bin->counts() / total_ontime;
537  bin->counts(rate);
538  }
539 
540  // Set background cube values from event cube
541  m_cube = eventcube.counts();
542 
543  } // endif: ontime was positive
544 
545  // Return
546  return;
547 
548 }
549 
550 
551 /***********************************************************************//**
552  * @brief Return spatially integrated background rate in units of
553  * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
554  *
555  * @param[in] logE Logarithm (base 10) of energy in TeV
556  * @return Spatially integrated background rate
557  * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
558  *
559  * Spatially integrates the background cube at a given energy. This method
560  * performs an interpolation between the energy maps. The integration is
561  * performed by summing over all bin contents multiplied by the solid angle
562  * of the bin.
563  ***************************************************************************/
564 double GCTACubeBackground::integral(const double& logE) const
565 {
566  // Update interpolation cache
567  update(logE);
568 
569  // Initialise result
570  double result = 0.0;
571 
572  // Loop over all map pixels
573  for (int i = 0; i < m_cube.npix(); ++i) {
574 
575  // Get bin value
576  double value = m_wgt_left * m_cube(i, m_inx_left) +
578 
579  // Sum bin contents
580  result += value * m_cube.solidangle(i);
581 
582  }
583 
584  // Return result
585  return result;
586 }
587 
588 
589 /***********************************************************************//**
590  * @brief Read background cube from FITS object
591  *
592  * @param[in] fits FITS object.
593  *
594  * Read the background cube from a FITS object.
595  ***************************************************************************/
597 {
598  // Clear object
599  clear();
600 
601  // Get HDUs
602  const GFitsImage& hdu_bgdcube = *fits.image("Primary");
603  const GFitsTable& hdu_ebounds = *fits.table(gammalib::extname_ebounds);
604 
605  // Read cube
606  m_cube.read(hdu_bgdcube);
607 
608  // Read energy boundaries
609  m_ebounds.read(hdu_ebounds);
610 
611  // Set energy node array
612  set_eng_axis();
613 
614  // Return
615  return;
616 }
617 
618 
619 /***********************************************************************//**
620  * @brief Write CTA background cube into FITS object.
621  *
622  * @param[in] fits FITS file.
623  ***************************************************************************/
625 {
626  // Write cube
627  m_cube.write(fits);
628 
629  // Write energy boundaries
630  m_ebounds.write(fits);
631 
632  // Return
633  return;
634 }
635 
636 
637 /***********************************************************************//**
638  * @brief Load background cube from FITS file
639  *
640  * @param[in] filename FITS file.
641  *
642  * Loads the background cube from a FITS file.
643  ***************************************************************************/
644 void GCTACubeBackground::load(const GFilename& filename)
645 {
646  // Put into OpenMP criticial zone
647  #pragma omp critical(GCTACubeBackground_load)
648  {
649 
650  // Open FITS file
651  GFits fits(filename);
652 
653  // Read background cube
654  read(fits);
655 
656  // Close FITS file
657  fits.close();
658 
659  } // end of OpenMP critical zone
660 
661  // Store filename
663 
664  // Return
665  return;
666 
667 }
668 
669 
670 /***********************************************************************//**
671  * @brief Save background cube into FITS file
672  *
673  * @param[in] filename background cube FITS file name.
674  * @param[in] clobber Overwrite existing file?
675  *
676  * Save the background cube into a FITS file.
677  ***************************************************************************/
678 void GCTACubeBackground::save(const GFilename& filename,
679  const bool& clobber) const
680 {
681  // Put into OpenMP criticial zone
682  #pragma omp critical(GCTACubeBackground_save)
683  {
684 
685  // Open or create FITS file
686  GFits fits;
687 
688  // Write background cube
689  write(fits);
690 
691  // Save FITS file
692  fits.saveto(filename, clobber);
693 
694  // Close Edisp file
695  fits.close();
696 
697  } // end of OpenMP critical zone
698 
699  // Store filename
701 
702  // Return
703  return;
704 }
705 
706 
707 /***********************************************************************//**
708  * @brief Print background information
709  *
710  * @param[in] chatter Chattiness.
711  * @return String containing background information.
712  ***************************************************************************/
713 std::string GCTACubeBackground::print(const GChatter& chatter) const
714 {
715  // Initialise result string
716  std::string result;
717 
718  // Continue only if chatter is not silent
719  if (chatter != SILENT) {
720 
721  // Append header
722  result.append("=== GCTACubeBackground ===");
723 
724  // Append parameters
725  result.append("\n"+gammalib::parformat("Filename")+m_filename);
726 
727  // Append energies
728  if (m_ebounds.size() > 0) {
729  result.append("\n"+m_ebounds.print(chatter));
730  }
731  else {
732  result.append("\n"+gammalib::parformat("Energy boundaries") +
733  "Not defined");
734  }
735 
736  // Append skymap definition
737  result.append("\n"+m_cube.print(chatter));
738 
739  } // endif: chatter was not silent
740 
741  // Return result
742  return result;
743 }
744 
745 
746 /*==========================================================================
747  = =
748  = Private methods =
749  = =
750  ==========================================================================*/
751 
752 /***********************************************************************//**
753  * @brief Initialise class members
754  ***************************************************************************/
756 {
757  // Initialise members
758  m_filename.clear();
759  m_cube.clear();
760  m_ebounds.clear();
761  m_elogmeans.clear();
762 
763  // Initialise cache
764  m_inx_left = 0;
765  m_inx_right = 0;
766  m_wgt_left = 0.0;
767  m_wgt_right = 0.0;
768 
769  // Return
770  return;
771 }
772 
773 
774 /***********************************************************************//**
775  * @brief Copy class members
776  *
777  * @param[in] bgd Background.
778  ***************************************************************************/
780 {
781  // Copy members
782  m_filename = bgd.m_filename;
783  m_cube = bgd.m_cube;
784  m_ebounds = bgd.m_ebounds;
785  m_elogmeans = bgd.m_elogmeans;
786 
787  // Copy cache
788  m_inx_left = bgd.m_inx_left;
789  m_inx_right = bgd.m_inx_right;
790  m_wgt_left = bgd.m_wgt_left;
791  m_wgt_right = bgd.m_wgt_right;
792 
793  // Return
794  return;
795 }
796 
797 
798 /***********************************************************************//**
799  * @brief Delete class members
800  ***************************************************************************/
802 {
803  // Return
804  return;
805 }
806 
807 
808 /***********************************************************************//**
809  * @brief Set nodes for a logarithmic (base 10) energy axis
810  *
811  *
812  * Set axis nodes so that each node is the logarithm of the energy values.
813  ***************************************************************************/
815 {
816  // Get number of bins
817  int bins = m_ebounds.size();
818 
819  // Clear node array
820  m_elogmeans.clear();
821 
822  // Compute nodes
823  for (int i = 0; i < bins; ++i) {
824 
825  // Get logE/TeV
827 
828  } // endfor: looped over energy bins
829 
830  // Return
831  return;
832 }
833 
834 
835 /***********************************************************************//**
836  * @brief Update 1D cache
837  *
838  * @param[in] logE Log10 energy in TeV.
839  *
840  * Updates the 1D interpolation cache. The interpolation cache is composed
841  * of two indices and weights that define 2 data values of the 2D skymap
842  * that are used for linear interpolation.
843  *
844  * @todo Write down formula
845  ***************************************************************************/
846 void GCTACubeBackground::update(const double& logE) const
847 {
848  // Set value for node array
849  m_elogmeans.set_value(logE);
850 
851  // Set indices and weighting factors for interpolation
856 
857  // Return
858  return;
859 }
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:472
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:1163
const GFilename & filename(void) const
Return background cube filename.
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
bool overlaps(const GSkyRegion &region) const
Checks whether a region overlaps with this map.
Definition: GSkyMap.cpp:1944
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
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:1192
const double & radius(void) const
Returns radius of region of interest in degrees.
Definition: GCTARoi.hpp:124
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
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.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Interface for the CTA region of interest class.
Definition: GCTARoi.hpp:48
void read(const GFitsTable &table)
Read energy boundaries from FITS table.
Definition: GEbounds.cpp:776
Time class.
Definition: GTime.hpp:54
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:580
void pointing(const GCTAPointing &pointing)
Set CTA pointing.
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2434
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:831
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:1267
std::string print(const GChatter &chatter=NORMAL) const
Print sky map.
Definition: GSkyMap.cpp:2557
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:275
Energy boundaries container class.
Definition: GEbounds.hpp:60
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:378
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:2467
Filename class.
Definition: GFilename.hpp:62
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
GGti gti(void) const
Get Good Time Intervals.
const GCTAInstDir & centre(void) const
Returns region of interest centre.
Definition: GCTARoi.hpp:110
virtual std::string instrument(void) const
Return instrument name.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
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:214
std::string print(const GChatter &chatter=NORMAL) const
Print energy boundaries.
Definition: GEbounds.cpp:1244
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:1140
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1090
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:269
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:245
void name(const std::string &name)
Set observation name.
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:360
double plaw_integral(const double &x1, const double &f1, const double &x2, const double &f2)
Returns the integral of a power law.
Definition: GMath.cpp:566
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.
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:151
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:1142
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
CTA observation class.
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
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:334
Observations container class interface definition.
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
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.