GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTABackground3D.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTABackground3D.cpp - CTA 3D background class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2014-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 GCTABackground3D.cpp
23  * @brief CTA 3D background 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 "GException.hpp"
33 #include "GFilename.hpp"
34 #include "GMath.hpp"
35 #include "GRan.hpp"
36 #include "GFits.hpp"
37 #include "GFitsBinTable.hpp"
38 #include "GCTAInstDir.hpp"
39 #include "GCTABackground3D.hpp"
40 #include "GCTASupport.hpp"
41 
42 /* __ Method name definitions ____________________________________________ */
43 #define G_MC "GCTABackground3D::mc(GEnergy&, GTime&, GRan&)"
44 #define G_SET_MEMBERS "GCTABackground3D::set_members()"
45 #define G_INIT_MC_CACHE "GCTABackground3D::init_mc_cache()"
46 
47 /* __ Macros _____________________________________________________________ */
48 
49 /* __ Coding definitions _________________________________________________ */
50 
51 /* __ Debug definitions __________________________________________________ */
52 //#define G_DEBUG_MC_INIT
53 //#define G_DEBUG_CACHE
54 
55 /* __ Constants __________________________________________________________ */
56 
57 
58 /*==========================================================================
59  = =
60  = Constructors/destructors =
61  = =
62  ==========================================================================*/
63 
64 /***********************************************************************//**
65  * @brief Void constructor
66  *
67  * Constructs empty background.
68  ***************************************************************************/
70 {
71  // Initialise class members
72  init_members();
73 
74  // Return
75  return;
76 }
77 
78 
79 /***********************************************************************//**
80  * @brief File constructor
81  *
82  * @param[in] filename FITS file name.
83  *
84  * Constructs background from a FITS file.
85  ***************************************************************************/
88 {
89  // Initialise class members
90  init_members();
91 
92  // Load background from file
93  load(filename);
94 
95  // Return
96  return;
97 }
98 
99 
100 /***********************************************************************//**
101  * @brief Copy constructor
102  *
103  * @param[in] bgd Background.
104  *
105  * Constructs background by copying from another background.
106  ***************************************************************************/
108  GCTABackground(bgd)
109 {
110  // Initialise class members
111  init_members();
112 
113  // Copy members
114  copy_members(bgd);
115 
116  // Return
117  return;
118 }
119 
120 
121 /***********************************************************************//**
122  * @brief Destructor
123  *
124  * Destructs background.
125  ***************************************************************************/
127 {
128  // Free members
129  free_members();
130 
131  // Return
132  return;
133 }
134 
135 
136 /*==========================================================================
137  = =
138  = Operators =
139  = =
140  ==========================================================================*/
141 
142 /***********************************************************************//**
143  * @brief Assignment operator
144  *
145  * @param[in] bgd Background.
146  * @return Background.
147  *
148  * Assigns background.
149  ***************************************************************************/
151 {
152  // Execute only if object is not identical
153  if (this != &bgd) {
154 
155  // Copy base class members
156  this->GCTABackground::operator=(bgd);
157 
158  // Free members
159  free_members();
160 
161  // Initialise private members
162  init_members();
163 
164  // Copy members
165  copy_members(bgd);
166 
167  } // endif: object was not identical
168 
169  // Return this object
170  return *this;
171 }
172 
173 
174 /***********************************************************************//**
175  * @brief Return background rate in units of events MeV\f$^{-1}\f$
176  * s\f$^{-1}\f$ sr\f$^{-1}\f$
177  *
178  * @param[in] logE Log10 of the true photon energy (TeV).
179  * @param[in] detx Tangential X coordinate in nominal system (radians).
180  * @param[in] dety Tangential Y coordinate in nominal system (radians).
181  * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
182  *
183  * Returns the background rate for a given energy and detector coordinates.
184  * The rate is given per ontime. The operator interpolates linearly in
185  * ``DETX`` and ``DETY``, and logarithmically in energy.
186  *
187  * The operator assures that the background rate never becomes negative. For
188  * invalid background models or detector coordinates outside the range
189  * covered by the response cube, the operator returns a rate of zero.
190  ***************************************************************************/
191 double GCTABackground3D::operator()(const double& logE,
192  const double& detx,
193  const double& dety) const
194 {
195  // Initialise background rate
196  double rate = 0.0;
197 
198  // Continue only if background model is valid and if the DETX and DETY
199  // values are in validity range
200  if (is_valid() && (detx >= m_detx_min) && (detx <= m_detx_max) &&
201  (dety >= m_dety_min) && (dety <= m_dety_max)) {
202 
203  // Retrieve references to energy node array
204  const GNodeArray& energy_nodes = m_background.axis_nodes(m_inx_energy);
205 
206  // Set energy for node array interpolation
207  energy_nodes.set_value(logE);
208 
209  // Bi-linearly interpolate the rates in both energy layers
210  double rate_emin = this->rate(energy_nodes.inx_left(), detx, dety);
211  double rate_emax = this->rate(energy_nodes.inx_right(), detx, dety);
212 
213  // If both rates are positive then perform a logarithmic interpolation
214  // in energy
215  if (rate_emin > 0.0 && rate_emax > 0.0) {
216 
217  // Perform energy interpolation
218  rate = std::exp(energy_nodes.wgt_left() * std::log(rate_emin) +
219  energy_nodes.wgt_right() * std::log(rate_emax));
220 
221  // Make sure that background rate is not negative
222  if (rate < 0.0) {
223  rate = 0.0;
224  }
225 
226  } // endif: both rates were positive
227 
228  } // endif: background model was valid
229 
230  // Return background rate
231  return rate;
232 }
233 
234 
235 /*==========================================================================
236  = =
237  = Public methods =
238  = =
239  ==========================================================================*/
240 
241 /***********************************************************************//**
242  * @brief Clear background
243  *
244  * Clears background.
245  ***************************************************************************/
247 {
248  // Free class members
249  free_members();
251 
252  // Initialise members
254  init_members();
255 
256  // Return
257  return;
258 }
259 
260 
261 /***********************************************************************//**
262  * @brief Clone background
263  *
264  * @return Deep copy of background.
265  *
266  * Returns a pointer to a deep copy of the background.
267  ***************************************************************************/
269 {
270  return new GCTABackground3D(*this);
271 }
272 
273 
274 /***********************************************************************//**
275  * @brief Read background from FITS table
276  *
277  * @param[in] table FITS table.
278  *
279  * Reads the background form the FITS @p table. The following column names
280  * are mandatory:
281  *
282  * DETX_LO - DETX lower bin boundaries
283  * DETX_HI - DETX upper bin boundaries
284  * DETY_LO - DETY lower bin boundaries
285  * DETY_HI - DETY upper bin boundaries
286  * ENERG_LO - Energy lower bin boundaries
287  * ENERG_HI - Energy upper bin boundaries
288  * BKG - Background template (or BGD as legacy column name)
289  *
290  * The data are stored in the m_background member. The ``DETX`` and ``DETY``
291  * axes will be set to radians, the energy axis will be set to log10. The
292  * data are expected to be given in unit of events per MeV, second and
293  * steradian, where the time is the real ontime.
294  *
295  * This method ignores all column names that are not the mandatory column
296  * names in the FITS @p table.
297  ***************************************************************************/
299 {
300  // Clear response table
302 
303  // Copy response table and skip all columns that are not mandatory. This
304  // kludge is needed to get rid of additional columns in the HESS
305  // background model.
306  GFitsTable *ptr = table.clone(); // Make sure that table is of same type
307  for (int i = 0; i < table.ncols(); ++i) {
308  std::string colname = table[i]->name();
309  if ((colname != "DETX_LO") && (colname != "DETX_HI") &&
310  (colname != "DETY_LO") && (colname != "DETY_HI") &&
311  (colname != "ENERG_LO") && (colname != "ENERG_HI") &&
312  (colname != "BKG") && (colname != "BGD")) {
313  ptr->remove(colname);
314  }
315  }
316 
317  // Read background table from reduced FITS table. The background table
318  // is expected to provide the number of events per MeV, second and
319  // steradian, where the time is the real ontime.
320  m_background.read(*ptr);
321 
322  // Delete reduced FITS table
323  delete ptr;
324 
325  // Set class members
326  set_members();
327 
328  // Return
329  return;
330 }
331 
332 /***********************************************************************//**
333  * @brief Write background into FITS table
334  *
335  * @param[in] table FITS binary table.
336  *
337  * Writes background into a FITS binary @p table.
338  *
339  * @todo Add necessary keywords.
340  ***************************************************************************/
342 {
343  // Write background table
344  m_background.write(table);
345 
346  // Return
347  return;
348 }
349 
350 
351 /***********************************************************************//**
352  * @brief Load background from FITS file
353  *
354  * @param[in] filename FITS file name.
355  *
356  * Loads the background from a FITS file.
357  *
358  * If no extension name is given the method scans the `HDUCLASS` keywords
359  * of all extensions and loads the background from the first extension
360  * for which `HDUCLAS4=BKG_3D`.
361  *
362  * Otherwise, the background will be loaded from the `BACKGROUND`
363  * extension.
364  ***************************************************************************/
365 void GCTABackground3D::load(const GFilename& filename)
366 {
367  // Open FITS file
368  GFits fits(filename);
369 
370  // Get the default extension name. If no GADF compliant name was found
371  // then set the default extension name to "BACKGROUND".
372  std::string extname = gammalib::gadf_hduclas4(fits, "BKG_3D");
373  if (extname.empty()) {
375  }
376 
377  // Get background table
378  const GFitsTable& table = *fits.table(filename.extname(extname));
379 
380  // Read effective area from table
381  read(table);
382 
383  // Close FITS file
384  fits.close();
385 
386  // Store filename
388 
389  // Return
390  return;
391 }
392 
393 
394 /***********************************************************************//**
395  * @brief Save background into FITS file
396  *
397  * @param[in] filename Background table FITS file name.
398  * @param[in] clobber Overwrite existing file? (default: false)
399  *
400  * Saves background into a FITS file. If a file with the given @p filename
401  * does not yet exist it will be created, otherwise the method opens the
402  * existing file. The method will create a (or replace an existing)
403  * background extension. The extension name can be specified as part
404  * of the @p filename, or if no extension name is given, is assumed to be
405  * `BACKGROUND`.
406  *
407  * An existing file will only be modified if the @p clobber flag is set to
408  * true.
409  ***************************************************************************/
410 void GCTABackground3D::save(const GFilename& filename, const bool& clobber) const
411 {
412  // Get extension name
413  std::string extname = filename.extname(gammalib::extname_cta_background3d);
414 
415  // Open or create FITS file (without extension name since the requested
416  // extension may not yet exist in the file)
417  GFits fits(filename.url(), true);
418 
419  // Remove extension if it exists already
420  if (fits.contains(extname)) {
421  fits.remove(extname);
422  }
423 
424  // Create binary table
426 
427  // Write the background table
428  write(table);
429 
430  // Set binary table extension name
431  table.extname(extname);
432 
433  // Append table to FITS file
434  fits.append(table);
435 
436  // Save to file
437  fits.save(clobber);
438 
439  // Return
440  return;
441 }
442 
443 
444 /***********************************************************************//**
445  * @brief Returns MC instrument direction
446  *
447  * @param[in] energy Event energy.
448  * @param[in] time Event trigger time.
449  * @param[in,out] ran Random number generator.
450  * @return CTA instrument direction.
451  *
452  * Returns a Monte Carlo simulated instrument direction for a given @p energy
453  * and event trigger @p time. The simulation is done using a rejection
454  * method that makes use of the background rate access operator. This assures
455  * that the simulation is consistent with the background rate interpolation
456  * that is done by the access operator.
457  ***************************************************************************/
459  const GTime& time,
460  GRan& ran) const
461 {
462  // If Monte Carlo has not been initialised then initialise the cache now
463  if (m_mc_max.empty()) {
464  init_mc_cache();
465  }
466 
467  // Determine energy range of response table
468  GEnergy emin = m_mc_spectrum.energy(0);
470  if (energy < emin || energy > emax) {
471  std::string msg = "Event energy "+energy.print()+" is outside the "
472  "energy range ["+emin.print()+", "+emax.print()+"] "
473  "covered by the background response table. Please "
474  "restrict the energy range of the simulation to the "
475  "validity range of the background response table.";
476  throw GException::invalid_value(G_MC, msg);
477  }
478 
479  // Allocate instrument direction
480  GCTAInstDir dir;
481 
482  // Continue only if there is a background model
483  if (!m_mc_max.empty()) {
484 
485  // Get reference to node array for the energy axis
486  const GNodeArray& energy_nodes = m_background.axis_nodes(m_inx_energy);
487 
488  // Get log10(energy) in TeV
489  double logE = energy.log10TeV();
490 
491  // Set values for node arrays
492  energy_nodes.set_value(logE);
493 
494  // Get indices for energy interpolation
495  int inx_left = energy_nodes.inx_left();
496  int inx_right = energy_nodes.inx_right();
497 
498  // Get maximum background rate as the maximum of the left and the right
499  // energy node. Add some margin.
500  double max_rate = m_mc_max[inx_left];
501  if (m_mc_max[inx_right] > max_rate) {
502  max_rate = m_mc_max[inx_right];
503  }
504  max_rate *= 1.5;
505 
506  // Get detx and dety widths
507  double detx_width = m_detx_max - m_detx_min;
508  double dety_width = m_dety_max - m_dety_min;
509 
510  // Get instrument direction using a rejection method. This assures
511  // that the Monte Carlo sample follows the model distribution
512  while (true) {
513 
514  // Get randomized detx and dety in radians
515  double detx = ran.uniform() * detx_width + m_detx_min;
516  double dety = ran.uniform() * dety_width + m_dety_min;
517 
518  // Get background rate for these coordinates. If the background
519  // rate is zero then get other coordinates.
520  double value = this->operator()(logE, detx, dety);
521  if (value <= 0.0) {
522  continue;
523  }
524 
525  // Dump a warning if value is larger than the maximum rate. This
526  // should never happen!
527  if (value > max_rate) {
528  std::string msg = "Background rate "+gammalib::str(value)+" "
529  "for "+energy.print()+" at "
530  "DETXY=["+
532  gammalib::str(dety*gammalib::rad2deg)+"] "
533  "is larger than the maximum expected rate "+
534  gammalib::str(max_rate)+". Something is "
535  "wrong with the code.";
536  gammalib::warning(G_MC, msg);
537  }
538 
539  // Get uniform random number
540  double uniform = ran.uniform() * max_rate;
541 
542  // Exit loop if we're not larger than the background rate value
543  if (uniform <= value) {
544  dir.detx(detx);
545  dir.dety(dety);
546  break;
547  }
548 
549  } // endwhile: loop until instrument direction was accepted
550 
551  } // endif: there were cube pixels and maps
552 
553  // Return instrument direction
554  return dir;
555 }
556 
557 
558 /***********************************************************************//**
559  * @brief Returns background rate integrated over energy interval in units
560  * of events s\f$^{-1}\f$ sr\f$^{-1}\f$
561  *
562  * @param[in] dir Instrument direction.
563  * @param[in] emin Minimum energy of energy interval.
564  * @param[in] emax Maximum energy of energy interval.
565  * @return Integrated background count rate (events s\f$^{-1}\f$ sr\f$^{-1}\f$).
566  *
567  * Returns the background count rate for a given instrument direction that
568  * is integrated over a specified energy interval. The rate is given per
569  * ontime.
570  *
571  * If the energy interval is not positive, a zero background rate is
572  * returned.
573  ***************************************************************************/
575  const GEnergy& emin,
576  const GEnergy& emax) const
577 {
578  // Initialise rate
579  double rate = 0.0;
580 
581  // Continue only if the background template is valid, the DETX and DETY
582  // values are comprised in template range and the energy interval is
583  // positive
584  if (is_valid() && (dir.detx() >= m_detx_min) &&
585  (dir.detx() <= m_detx_max) &&
586  (dir.dety() >= m_dety_min) &&
587  (dir.dety() <= m_dety_max) &&
588  (emax > emin)) {
589 
590  // Initialise first and second node
591  double x1 = emin.MeV();
592  double x2 = 0.0;
593  double f1 = (*this)(emin.log10TeV(), dir.detx(), dir.dety());
594  double f2 = 0.0;
595 
596  // Loop over all nodes
597  for (int i = 0; i < m_energy.size(); ++i) {
598 
599  // If node energy is below x1 then skip node
600  if (m_energy[i].MeV() <= x1) {
601  continue;
602  }
603 
604  // If node energy is above emax then use emax as energy
605  if (m_energy[i] > emax) {
606  x2 = emax.MeV();
607  f2 = (*this)(emax.log10TeV(), dir.detx(), dir.dety());
608  }
609 
610  // ... otherwise use node energy
611  else {
612  x2 = m_energy[i].MeV();
613  f2 = this->rate(i, dir.detx(), dir.dety());
614  }
615 
616  // Compute integral
617  if (f1 > 0.0 && f2 > 0.0) {
618  rate += gammalib::plaw_integral(x1, f1, x2, f2);
619  }
620 
621  // Set second node as first node
622  x1 = x2;
623  f1 = f2;
624 
625  // If node energy is above emax then break now
626  if (m_energy[i] > emax) {
627  break;
628  }
629 
630  } // endfor: looped over all nodes
631 
632  // If last node energy is below emax then compute last part of
633  // integral up to emax
634  if (x1 < emax.MeV()) {
635  x2 = emax.MeV();
636  f2 = (*this)(emax.log10TeV(), dir.detx(), dir.dety());
637  if (f1 > 0.0 && f2 > 0.0) {
638  rate += gammalib::plaw_integral(x1, f1, x2, f2);
639  }
640  }
641 
642  } // endif: energy interval was positive
643 
644  // Return background rate
645  return rate;
646 }
647 
648 
649 /***********************************************************************//**
650  * @brief Print background information
651  *
652  * @param[in] chatter Chattiness.
653  * @return String containing background information.
654  ***************************************************************************/
655 std::string GCTABackground3D::print(const GChatter& chatter) const
656 {
657  // Initialise result string
658  std::string result;
659 
660  // Continue only if chatter is not silent
661  if (chatter != SILENT) {
662 
663  // Append header
664  result.append("=== GCTABackground3D ===");
665 
666  // Append information
667  result.append("\n"+gammalib::parformat("Filename")+m_filename);
668 
669  // If background model is valid then print information
670  if (is_valid()) {
671 
672  // Compute DETX boundaries in deg
673  double detx_min = m_background.axis_lo(m_inx_detx, 0);
674  double detx_max = m_background.axis_hi(m_inx_detx, m_num_detx-1);
675 
676  // Compute DETY boundaries in deg
677  double dety_min = m_background.axis_lo(m_inx_dety,0);
678  double dety_max = m_background.axis_hi(m_inx_dety, m_num_dety-1);
679 
680  // Compute energy boundaries in TeV
681  double emin = m_background.axis_lo(m_inx_energy, 0);
682  double emax = m_background.axis_hi(m_inx_energy, m_num_energy-1);
683 
684  // Append information
685  result.append("\n"+gammalib::parformat("Number of DETX bins") +
687  result.append("\n"+gammalib::parformat("Number of DETY bins") +
689  result.append("\n"+gammalib::parformat("Number of energy bins") +
691  result.append("\n"+gammalib::parformat("DETX range"));
692  result.append(gammalib::str(detx_min)+" - "+gammalib::str(detx_max)+" deg");
693  result.append("\n"+gammalib::parformat("DETX range"));
694  result.append(gammalib::str(dety_min)+" - "+gammalib::str(dety_max)+" deg");
695  result.append("\n"+gammalib::parformat("Energy range"));
696  result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
697 
698  } // endif: there were 3 axis
699 
700  // ... otherwise show empty array
701  else {
702  result.append("\n"+gammalib::parformat("Number of DETX bins") +
703  gammalib::str(0));
704  result.append("\n"+gammalib::parformat("Number of DETY bins") +
705  gammalib::str(0));
706  result.append("\n"+gammalib::parformat("Number of energy bins") +
707  gammalib::str(0));
708  }
709 
710  } // endif: chatter was not silent
711 
712  // Return result
713  return result;
714 }
715 
716 
717 /*==========================================================================
718  = =
719  = Private methods =
720  = =
721  ==========================================================================*/
722 
723 /***********************************************************************//**
724  * @brief Initialise class members
725  ***************************************************************************/
727 {
728  // Initialise members
729  m_filename.clear();
731  m_energy.clear();
732  m_inx_detx = 0;
733  m_inx_dety = 1;
734  m_inx_energy = 2;
735  m_inx_bgd = 0;
736  m_num_detx = 0;
737  m_num_dety = 0;
738  m_num_energy = 0;
739  m_num[0] = 0;
740  m_num[1] = 0;
741  m_num[2] = 0;
742  m_detx_min = 0.0;
743  m_detx_max = 0.0;
744  m_dety_min = 0.0;
745  m_dety_max = 0.0;
746  m_logE_min = 0.0;
747  m_logE_max = 0.0;
748 
749  // Initialise MC cache
750  m_mc_max.clear();
752 
753  // Return
754  return;
755 }
756 
757 
758 /***********************************************************************//**
759  * @brief Copy class members
760  *
761  * @param[in] bgd Background.
762  ***************************************************************************/
764 {
765  // Copy members
766  m_filename = bgd.m_filename;
768  m_energy = bgd.m_energy;
769  m_inx_detx = bgd.m_inx_detx;
770  m_inx_dety = bgd.m_inx_dety;
772  m_inx_bgd = bgd.m_inx_bgd;
773  m_num_detx = bgd.m_num_detx;
774  m_num_dety = bgd.m_num_dety;
776  m_num[0] = bgd.m_num[0];
777  m_num[1] = bgd.m_num[1];
778  m_num[2] = bgd.m_num[2];
779  m_detx_min = bgd.m_detx_min;
780  m_detx_max = bgd.m_detx_max;
781  m_dety_min = bgd.m_dety_min;
782  m_dety_max = bgd.m_dety_max;
783  m_logE_min = bgd.m_logE_min;
784  m_logE_max = bgd.m_logE_max;
785 
786  // Copy MC cache
787  m_mc_max = bgd.m_mc_max;
789 
790  // Return
791  return;
792 }
793 
794 
795 /***********************************************************************//**
796  * @brief Delete class members
797  ***************************************************************************/
799 {
800  // Return
801  return;
802 }
803 
804 
805 /***********************************************************************//**
806  * @brief Set members from background table
807  *
808  * @exception GException::invalid_value
809  * Response table is not three-dimensional.
810  *
811  * Set class members based on the background table. The DETX and DETY axes
812  * will be set to radians, the energy axis will be set to log10.
813  ***************************************************************************/
815 {
816  // Throw an exception if the table is not three-dimensional
817  if (m_background.axes() != 3) {
818  std::string msg = "Expected three-dimensional background "
819  "response table but found "+
821  " dimensions. Please specify a three-dimensional "
822  "background.";
824  }
825 
826  // Get mandatory indices (throw exception if not found)
827  m_inx_detx = m_background.axis("DETX");
828  m_inx_dety = m_background.axis("DETY");
829  m_inx_energy = m_background.axis("ENERG");
830  if (m_background.has_table("BKG")) {
831  m_inx_bgd = m_background.table("BKG");
832  }
833  else {
834  m_inx_bgd = m_background.table("BGD"); // Old name, should not be used
835  }
836 
837  // Get axes dimensions
841 
842  // Set dimension array
846 
847  // Set DETX and DETY axis to radians
850 
851  // Set energy axis to logarithmic scale
853 
854  // Compute DETX boundaries in radians
859 
860  // Compute DETY boundaries in radians
865 
866  // Compute energy boundaries in log10(TeV)
871  m_logE_min = emin.log10TeV();
872  m_logE_max = emax.log10TeV();
873 
874  // Setup energy vector
875  m_energy.clear();
876  for (int i = 0; i < m_num_energy; ++i) {
881  GEnergy energy(std::sqrt(elo.TeV()*ehi.TeV()), "TeV");
882  m_energy.append(energy);
883  }
884 
885  // Return
886  return;
887 }
888 
889 
890 /***********************************************************************//**
891  * @brief Return background rate bin index
892  *
893  * @param[in] idetx DETX index.
894  * @param[in] idety DETY index.
895  * @param[in] iebin Energy index.
896  * @return Background rate bin index.
897  *
898  * Returns the background rate bin index independent of the ordering.
899  ***************************************************************************/
900 int GCTABackground3D::index(const int& idetx,
901  const int& idety,
902  const int& iebin) const
903 {
904  // Set index arrays
905  int inx[3];
906  inx[m_inx_detx] = idetx;
907  inx[m_inx_dety] = idety;
908  inx[m_inx_energy] = iebin;
909 
910  // Compute index
911  int index = inx[0] + (inx[1] + inx[2] * m_num[1]) * m_num[0];
912 
913  // Return index
914  return index;
915 }
916 
917 
918 /***********************************************************************//**
919  * @brief Initialise Monte Carlo cache
920  *
921  * @exception GException::invalid_value
922  * Background response table was empty.
923  *
924  * Initialises the cache for Monte Carlo sampling.
925  ***************************************************************************/
927 {
928  // Set number of subbins per energy bin of response cube
929  const int nsubbins = 10;
930 
931  // Initialise cache
932  m_mc_max.clear();
934 
935  // Initialise maximum rate array
937 
938  // Compute DETX and DETY binsize in radians
939  double detx_bin = (m_detx_max - m_detx_min) / m_num_detx;
940  double dety_bin = (m_dety_max - m_dety_min) / m_num_dety;
941 
942  // Set energy bins for Monte Carlo cache spectrum. Make sure that the
943  // actual energy nodes are part of the energy bins to avoid that sharp
944  // structures in the background rate variations are missed.
945  GEnergies energies;
946  for (int i = 0; i < m_num_energy; ++i) {
947 
948  // If this is the first node then add energies below the first node
949  // down to the lower energy limit of the background response
950  if (i == 0) {
953  GEnergies enodes(nsubbins/2, emin, m_energy[i]);
954  enodes.remove(enodes.size()-1); // avoids duplication of energy
955  energies.extend(enodes);
956  }
957 
958  // If this is not the last node then append nodes between current
959  // node and next node
960  if (i < m_num_energy-1) {
961  GEnergies enodes(nsubbins, m_energy[i], m_energy[i+1]);
962  enodes.remove(enodes.size()-1); // avoids duplication of energy
963  energies.extend(enodes);
964  }
965 
966  // ... otherwise this is the last node. Add energies above the last
967  // node up to the upper energy limit of the background response
968  else {
971  GEnergies enodes(nsubbins/2, m_energy[i], emax);
972  energies.extend(enodes);
973  }
974 
975  } // endfor: looped over energy bins of background response
976 
977  // Loop over all energy bins
978  for (int i = 0; i < energies.size(); ++i) {
979 
980  // Get logE of energy bin
981  double logE = energies[i].log10TeV();
982 
983  // Initialise cache with cumulative pixel fluxes and compute total flux
984  // in response table for normalization. Negative pixels are excluded
985  // from the cumulative map.
986  double total_flux = 0.0; // units: events/s/MeV
987  double dx = 0.5 * detx_bin;
988  double dy = 0.5 * dety_bin;
989  double detx = m_detx_min + dx;
990  for (int ix = 0; ix < m_num_detx; ++ix, detx += detx_bin) {
991  double dety = m_dety_min + dy;
992  for (int iy = 0; iy < m_num_dety; ++iy, dety += dety_bin) {
993 
994  // Compute intensities
995  double i0 = (*this)(logE, detx, dety);
996  double i1 = (*this)(logE, detx-dx, dety-dy);
997  double i2 = (*this)(logE, detx, dety-dy);
998  double i3 = (*this)(logE, detx+dx, dety-dy);
999  double i4 = (*this)(logE, detx+dx, dety);
1000  double i5 = (*this)(logE, detx+dx, dety+dy);
1001  double i6 = (*this)(logE, detx, dety+dy);
1002  double i7 = (*this)(logE, detx-dx, dety+dy);
1003  double i8 = (*this)(logE, detx-dx, dety);
1004 
1005  // Compute solid angle of the 8 pixel wedges
1006  double s1 = solid_angle(detx, dety, detx-dx, dety-dy, detx, dety-dy);
1007  double s2 = solid_angle(detx, dety, detx, dety-dy, detx+dx, dety-dy);
1008  double s3 = solid_angle(detx, dety, detx+dx, dety-dy, detx+dx, dety);
1009  double s4 = solid_angle(detx, dety, detx+dx, dety, detx+dx, dety+dy);
1010  double s5 = solid_angle(detx, dety, detx+dx, dety+dy, detx, dety+dy);
1011  double s6 = solid_angle(detx, dety, detx, dety+dy, detx-dx, dety+dy);
1012  double s7 = solid_angle(detx, dety, detx-dx, dety+dy, detx-dx, dety);
1013  double s8 = solid_angle(detx, dety, detx-dx, dety, detx-dx, dety-dy);
1014 
1015  // Compute flux by summing the flux in 8 pixel wedges
1016  double flux1 = s1 * (i1 + i2 + i0);
1017  double flux2 = s2 * (i2 + i3 + i0);
1018  double flux3 = s3 * (i3 + i4 + i0);
1019  double flux4 = s4 * (i4 + i5 + i0);
1020  double flux5 = s5 * (i5 + i6 + i0);
1021  double flux6 = s6 * (i6 + i7 + i0);
1022  double flux7 = s7 * (i7 + i8 + i0);
1023  double flux8 = s8 * (i8 + i1 + i0);
1024  double flux = (flux1 + flux2 + flux3 + flux4 +
1025  flux5 + flux6 + flux7 + flux8) / 3.0;
1026 
1027  // Sum flux
1028  if (flux > 0.0) {
1029  total_flux += flux;
1030  }
1031 
1032  } // endfor: loop over DETY
1033 
1034  } // endfor: loop over DETX
1035 
1036  // Append node if flux is positive
1037  if (total_flux > 0.0) {
1038  m_mc_spectrum.append(energies[i], total_flux);
1039  }
1040 
1041  // Dump spectrum for debugging
1042  #if defined(G_DEBUG_MC_INIT)
1043  std::cout << "Energy=" << energies[i];
1044  std::cout << " Rate_total=" << total_flux;
1045  std::cout << " events/(s MeV)" << std::endl;
1046  #endif
1047 
1048  } // endfor: looped over energy bins
1049 
1050  // Make sure that spectrum is not empty (if this is the case the entire
1051  // background response table was empty)
1052  if (m_mc_spectrum.nodes() == 0) {
1053  std::string msg = "Background response table is empty. Please provide "
1054  "a valid three-dimensional background template.";
1056  }
1057 
1058  // Return
1059  return;
1060 }
1061 
1062 
1063 /***********************************************************************//**
1064  * @brief Initialise array of maximum background rate
1065  *
1066  * Initialises the array m_mc_max that holds the maximum background rate in
1067  * each DETX-DETY plane of the background model.
1068  ***************************************************************************/
1070 {
1071  // Initialise cache
1072  m_mc_max.clear();
1073 
1074  // Loop over all maps to determine the maximum map intensity in units
1075  // of events/s/sr/MeV
1076  for (int iebin = 0; iebin < m_num_energy; ++iebin) {
1077 
1078  // Initialise maximum rate
1079  double max_rate = 0.0;
1080 
1081  // Loop over DETX-DETY plane
1082  for (int ix = 0; ix < m_num_detx; ++ix) {
1083  for (int iy = 0; iy < m_num_dety; ++iy) {
1084 
1085  // Get bin index
1086  int inx = index(ix, iy, iebin);
1087 
1088  // Get background rate
1089  double rate = m_background(m_inx_bgd, inx);
1090 
1091  // If background rate is larger than maximum then store the
1092  // background rate as new maximum
1093  if (rate > max_rate) {
1094  max_rate = rate;
1095  }
1096 
1097  } // endfor: looped over DETY
1098  } // endfor: looped over DETX
1099 
1100  // Append maximum rate
1101  m_mc_max.push_back(max_rate);
1102 
1103  } // endif: looped over all maps
1104 
1105  // Return
1106  return;
1107 }
1108 
1109 
1110 /***********************************************************************//**
1111  * @brief Compute solid angle of pixel wedge
1112  *
1113  * @param[in] detx1 DETX of first edge.
1114  * @param[in] dety1 DETY of first edge.
1115  * @param[in] detx2 DETX of second edge.
1116  * @param[in] dety2 DETY of second edge.
1117  * @param[in] detx3 DETX of third edge.
1118  * @param[in] dety3 DETY of third edge.
1119  * @return Solid angle (steradians).
1120  *
1121  * Estimate the solid angle subtended by 3 coordinates using Huilier's
1122  * theorem.
1123  *
1124  * Below, the definiton of the pixel cornes and sides are shown as used
1125  * within the code.
1126  *
1127  * a12
1128  * 1---------2
1129  * | /
1130  * | /
1131  * | /
1132  * | /
1133  * a13| /a23
1134  * | /
1135  * | /
1136  * | /
1137  * |/
1138  * 3
1139  *
1140  ***************************************************************************/
1141 double GCTABackground3D::solid_angle(const double& detx1, const double& dety1,
1142  const double& detx2, const double& dety2,
1143  const double& detx3, const double& dety3) const
1144 {
1145  // Initialise solid angle
1146  double solidangle = 0.0;
1147 
1148  // Set sky directions
1149  double theta1 = std::sqrt(detx1 * detx1 + dety1 * dety1);
1150  double phi1 = std::atan2(dety1, detx1);
1151  double theta2 = std::sqrt(detx2 * detx2 + dety2 * dety2);
1152  double phi2 = std::atan2(dety2, detx2);
1153  double theta3 = std::sqrt(detx3 * detx3 + dety3 * dety3);
1154  double phi3 = std::atan2(dety3, detx3);
1155  GSkyDir dir1;
1156  GSkyDir dir2;
1157  GSkyDir dir3;
1158  dir1.radec(phi1, gammalib::pihalf - theta1);
1159  dir2.radec(phi2, gammalib::pihalf - theta2);
1160  dir3.radec(phi3, gammalib::pihalf - theta3);
1161 
1162  // Compute angular distances between pixel corners
1163  double a12 = dir1.dist(dir2);
1164  double a13 = dir1.dist(dir3);
1165  double a23 = dir2.dist(dir3);
1166 
1167  // Compute solid angle
1168  double s = 0.5 * (a12 + a23 + a13);
1169  solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
1170  std::tan(0.5*(s-a12)) *
1171  std::tan(0.5*(s-a23)) *
1172  std::tan(0.5*(s-a13))));
1173 
1174  // Return solid angle
1175  return solidangle;
1176 }
1177 
1178 
1179 /***********************************************************************//**
1180  * @brief Return background rate for a given energy bin and DETX-DETY value
1181  * (events/s/MeV/sr)
1182  *
1183  * @param[in] iebin Energy index.
1184  * @param[in] detx Tangential X coordinate in nominal system (radians).
1185  * @param[in] dety Tangential Y coordinate in nominal system (radians).
1186  * @return Background rate (events/s/MeV/sr)
1187  ***************************************************************************/
1188 double GCTABackground3D::rate(const int& iebin, const double& detx,
1189  const double& dety) const
1190 {
1191  // Retrieve references to node arrays
1192  const GNodeArray& detxs = m_background.axis_nodes(m_inx_detx);
1193  const GNodeArray& detys = m_background.axis_nodes(m_inx_dety);
1194 
1195  // Set DETX and DETY values for node array interpolation
1196  detxs.set_value(detx);
1197  detys.set_value(dety);
1198 
1199  // Get bin indices
1200  int inx_ll = index(detxs.inx_left(), detys.inx_left(), iebin);
1201  int inx_lr = index(detxs.inx_left(), detys.inx_right(), iebin);
1202  int inx_rl = index(detxs.inx_right(), detys.inx_left(), iebin);
1203  int inx_rr = index(detxs.inx_right(), detys.inx_right(), iebin);
1204 
1205  // Bi-linearly interpolate the rates
1206  double rate = detxs.wgt_left() *
1207  (m_background(m_inx_bgd, inx_ll) * detys.wgt_left() +
1208  m_background(m_inx_bgd, inx_lr) * detys.wgt_right()) +
1209  detxs.wgt_right() *
1210  (m_background(m_inx_bgd, inx_rl) * detys.wgt_left() +
1211  m_background(m_inx_bgd, inx_rr) * detys.wgt_right());
1212 
1213  // Make sure that background rate is not negative
1214  if (rate < 0.0) {
1215  rate = 0.0;
1216  }
1217 
1218  // Return
1219  return rate;
1220 }
double m_dety_max
DETY maximum (radians)
const GNodeArray & axis_nodes(const int &axis) const
Return axis nodes.
const int & ncols(void) const
Return number of columns in table.
Definition: GFitsTable.hpp:134
void clear(void)
Clear energy container.
Definition: GEnergies.cpp:225
GCTABackground3D & operator=(const GCTABackground3D &bgd)
Assignment operator.
GFilename m_filename
Name of background response file.
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
void init_mc_max_rate(void) const
Initialise array of maximum background rate.
Node array class.
Definition: GNodeArray.hpp:60
#define G_INIT_MC_CACHE
void detx(const double &x)
Set DETX coordinate (in radians)
double m_dety_min
DETY minimum (radians)
Random number generator class definition.
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition: GTools.cpp:1265
int m_num_dety
Number of DETY bins.
int m_num[3]
Array of number of bins.
const GCTAResponseTable & table(void) const
Return response table.
std::vector< double > m_mc_max
Maximum background rate.
GCTABackground & operator=(const GCTABackground &bgd)
Assignment operator.
void write(GFitsBinTable &table) const
Write background into FITS table.
const std::string extname_cta_background3d
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
GEnergies m_energy
Vector of energies.
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
Random number generator class.
Definition: GRan.hpp:44
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
double solid_angle(const double &detx1, const double &dety1, const double &detx2, const double &dety2, const double &detx3, const double &dety3) const
Compute solid angle of pixel wedge.
FITS file class.
Definition: GFits.hpp:63
void axis_radians(const int &axis)
Set nodes for a radians axis.
FITS file class interface definition.
CTA 3D background class.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
int m_inx_detx
DETX index.
void read(const GFitsTable &table)
Read response table from FITS table HDU.
Energy container class.
Definition: GEnergies.hpp:60
void clear(void)
Clear background.
const double & axis_lo(const int &axis, const int &bin) const
Return lower bin boundary for bin in axis.
void append(const GEnergy &energy, const double &intensity)
Append node.
Definition of support function used by CTA classes.
void write(GFitsTable &table) const
Write response table into FITS table HDU.
int axis(const std::string &name) const
Determine index of an axis.
void radec(const double &ra, const double &dec)
Set equatorial sky direction (radians)
Definition: GSkyDir.cpp:175
void copy_members(const GCTABackground3D &bgd)
Copy class members.
void remove(const int &index)
Remove energy from container.
Definition: GEnergies.cpp:359
const double pihalf
Definition: GMath.hpp:38
double m_detx_max
DETX maximum (radians)
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
#define G_SET_MEMBERS
int m_num_energy
Number of energy bins.
void init_members(void)
Initialise class members.
bool has_table(const std::string &name) const
Check whether a table exists.
void save(const GFilename &filename, const bool &clobber=false) const
Save background into FITS file.
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:275
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:848
void free_members(void)
Delete class members.
GVector tan(const GVector &vector)
Computes tangens of vector elements.
Definition: GVector.cpp:1289
int m_num_detx
Number of DETX bins.
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
#define G_MC
void set_members(void)
Set members from background table.
void clear(void)
Clear response table.
Filename class.
Definition: GFilename.hpp:62
GCTAInstDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Returns MC instrument direction.
int index(const int &idetx, const int &idety, const int &iebin) const
Return background rate bin index.
double m_logE_min
Log10(E/TeV) minimum.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:162
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void init_mc_cache(void) const
Initialise Monte Carlo cache.
CTA instrument direction class interface definition.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
bool is_valid(void) const
Return validity of background model.
GChatter
Definition: GTypemaps.hpp:33
GCTAResponseTable m_background
Background response table.
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
GEnergy energy(const int &index) const
Return node energy.
CTA 3D background class definition.
void extend(const GEnergies &energies)
Append energy container.
Definition: GEnergies.cpp:383
void load(const GFilename &filename)
Load background from FITS file.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
void axis_log10(const int &axis)
Set nodes for a logarithmic (base 10) axis.
int table(const std::string &name) const
Determine index of table.
std::string gadf_hduclas4(const GFits &fits, const std::string &hduclas4)
Return extension name for GADF response table of given HDU class 4.
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
void init_members(void)
Initialise class members.
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
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
int m_inx_dety
DETY index.
virtual double operator()(const double &logE, const double &detx, const double &dety) const
Return background rate in units of events MeV s sr .
int axis_bins(const int &axis) const
Return number bins in an axis.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
virtual void clear(void)
Clear spectral nodes model.
double m_logE_max
Log10(E/TeV) maximum.
GModelSpectralNodes m_mc_spectrum
Response cube spectrum.
int nodes(void) const
Return number of nodes.
int m_inx_energy
Energy index.
FITS binary table class.
const double & axis_hi(const int &axis, const int &bin) const
Return upper bin boundary for bin in axis.
void remove(const int &colnum)
Remove column from the table.
Definition: GFitsTable.cpp:516
virtual GFitsTable * clone(void) const =0
Clones object.
GCTABackground3D(void)
Void constructor.
GEnergy & append(const GEnergy &energy)
Append energy to container.
Definition: GEnergies.cpp:303
Exception handler interface definition.
Abstract base class for the CTA background model.
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:265
FITS binary table class definition.
double atan2(const double &y, const double &x)
Compute arc tangens in radians.
Definition: GMath.cpp:102
std::string print(const GChatter &chatter=NORMAL) const
Print background information.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
void dety(const double &y)
Set DETY coordinate (in radians)
const std::string & axis_lo_unit(const int &axis) const
Return lower bin boundary unit for axis.
double rate(const int &iebin, const double &detx, const double &dety) const
Return background rate for a given energy bin and DETX-DETY value (events/s/MeV/sr) ...
const std::string & axis_hi_unit(const int &axis) const
Return upper bin boundary unit for axis.
const double rad2deg
Definition: GMath.hpp:44
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
Sky direction class.
Definition: GSkyDir.hpp:62
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
const int & axes(void) const
Return number of axes of the tables.
double rate_ebin(const GCTAInstDir &dir, const GEnergy &emin, const GEnergy &emax) const
Returns background rate integrated over energy interval in units of events s sr .
GVector atan(const GVector &vector)
Computes arctan of vector elements.
Definition: GVector.cpp:1058
GFilename filename(void) const
Return filename.
Filename class interface definition.
virtual ~GCTABackground3D(void)
Destructor.
Mathematical function definitions.
int m_inx_bgd
Background index.
void read(const GFitsTable &table)
Read background from FITS table.
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
void free_members(void)
Delete class members.
GCTABackground3D * clone(void) const
Clone background.
double m_detx_min
DETX minimum (radians)