GammaLib 2.0.0
Loading...
Searching...
No Matches
GCTABackground2D.cpp
Go to the documentation of this file.
1/***************************************************************************
2 * GCTABackground2D.cpp - CTA 2D background class *
3 * ----------------------------------------------------------------------- *
4 * copyright (C) 2021 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 GCTABackground2D.cpp
23 * @brief CTA 2D 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 "GCTABackground2D.hpp"
40#include "GCTASupport.hpp"
41
42/* __ Method name definitions ____________________________________________ */
43#define G_MC "GCTABackground2D::mc(GEnergy&, GTime&, GRan&)"
44#define G_SET_MEMBERS "GCTABackground2D::set_members()"
45#define G_INIT_MC_CACHE "GCTABackground2D::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
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
91
92 // Load background from file
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 computes the offset angle
185 *
186 * \f[\theta = \sqrt{ {\rm DETX}^2 + {\rm DETY}^2 }\f]
187 *
188 * and interpolates linearly in \f$\theta\f$ and logarithmically in energy.
189 *
190 * The operator assures that the background rate never becomes negative. For
191 * invalid background models or detector coordinates outside the range
192 * covered by the response cube, the operator returns a rate of zero.
193 ***************************************************************************/
194double GCTABackground2D::operator()(const double& logE,
195 const double& detx,
196 const double& dety) const
197{
198 // Initialise background rate
199 double rate = 0.0;
200
201 // Continue only if background model is valid
202 if (is_valid()) {
203
204 // Compute theta in radians
205 double theta = std::sqrt(detx * detx + dety * dety);
206
207 // Continue only if theta is within validity range
208 if ((theta >= m_theta_min) && (theta <= m_theta_max)) {
209
210 // Retrieve references to energy node array
211 const GNodeArray& energy_nodes = m_background.axis_nodes(m_inx_energy);
212
213 // Set energy for node array interpolation
214 energy_nodes.set_value(logE);
215
216 // Bi-linearly interpolate the rates in both energy layers
217 double rate_emin = this->rate(energy_nodes.inx_left(), theta);
218 double rate_emax = this->rate(energy_nodes.inx_right(), theta);
219
220 // If both rates are positive then perform a logarithmic interpolation
221 // in energy
222 if (rate_emin > 0.0 && rate_emax > 0.0) {
223
224 // Perform energy interpolation
225 rate = std::exp(energy_nodes.wgt_left() * std::log(rate_emin) +
226 energy_nodes.wgt_right() * std::log(rate_emax));
227
228 // Make sure that background rate is not negative
229 if (rate < 0.0) {
230 rate = 0.0;
231 }
232
233 } // endif: both rates were positive
234
235 } // endif: theta was valid
236
237 } // endif: background model was valid
238
239 // Return background rate
240 return rate;
241}
242
243
244/*==========================================================================
245 = =
246 = Public methods =
247 = =
248 ==========================================================================*/
249
250/***********************************************************************//**
251 * @brief Clear background
252 *
253 * Clears background.
254 ***************************************************************************/
256{
257 // Free class members
258 free_members();
260
261 // Initialise members
263 init_members();
264
265 // Return
266 return;
267}
268
269
270/***********************************************************************//**
271 * @brief Clone background
272 *
273 * @return Deep copy of background.
274 *
275 * Returns a pointer to a deep copy of the background.
276 ***************************************************************************/
278{
279 return new GCTABackground2D(*this);
280}
281
282
283/***********************************************************************//**
284 * @brief Read background from FITS table
285 *
286 * @param[in] table FITS table.
287 *
288 * Reads the background form the FITS @p table. The following column names
289 * are mandatory:
290 *
291 * THETA_LO - THETA lower bin boundaries
292 * THETA_HI - THETA upper bin boundaries
293 * ENERG_LO - Energy lower bin boundaries
294 * ENERG_HI - Energy upper bin boundaries
295 * BKG - Background template
296 *
297 * The data are stored in the m_background member. The ``THETA`` axis will
298 * be set to radians, the energy axis will be set to log10. The data are
299 * expected to be given in unit of events per MeV, second and steradian,
300 * where the time is the real ontime.
301 *
302 * This method ignores all column names that are not the mandatory column
303 * names in the FITS @p table.
304 ***************************************************************************/
306{
307 // Clear response table
309
310 // Copy response table and skip all columns that are not mandatory
311 GFitsTable *ptr = table.clone(); // Make sure that table is of same type
312 for (int i = 0; i < table.ncols(); ++i) {
313 std::string colname = table[i]->name();
314 if ((colname != "THETA_LO") && (colname != "THETA_HI") &&
315 (colname != "ENERG_LO") && (colname != "ENERG_HI") &&
316 (colname != "BKG")) {
317 ptr->remove(colname);
318 }
319 }
320
321 // Read background table from reduced FITS table. The background table
322 // is expected to provide the number of events per MeV, second and
323 // steradian, where the time is the real ontime.
324 m_background.read(*ptr);
325
326 // Delete reduced FITS table
327 delete ptr;
328
329 // Set class members
330 set_members();
331
332 // Return
333 return;
334}
335
336/***********************************************************************//**
337 * @brief Write background into FITS table
338 *
339 * @param[in] table FITS binary table.
340 *
341 * Writes background into a FITS binary @p table.
342 *
343 * @todo Add necessary keywords.
344 ***************************************************************************/
346{
347 // Write background table
349
350 // Return
351 return;
352}
353
354
355/***********************************************************************//**
356 * @brief Load background from FITS file
357 *
358 * @param[in] filename FITS file name.
359 *
360 * Loads the background from a FITS file.
361 *
362 * If no extension name is given the method scans the `HDUCLASS` keywords
363 * of all extensions and loads the background from the first extension
364 * for which `HDUCLAS4=BKG_2D`.
365 *
366 * Otherwise, the background will be loaded from the `BKG` extension.
367 ***************************************************************************/
369{
370 // Open FITS file
371 GFits fits(filename);
372
373 // Get the default extension name. If no GADF compliant name was found
374 // then set the default extension name to "BKG".
375 std::string extname = gammalib::gadf_hduclas4(fits, "BKG_2D");
376 if (extname.empty()) {
378 }
379
380 // Get background table
381 const GFitsTable& table = *fits.table(filename.extname(extname));
382
383 // Read effective area from table
384 read(table);
385
386 // Close FITS file
387 fits.close();
388
389 // Store filename
391
392 // Return
393 return;
394}
395
396
397/***********************************************************************//**
398 * @brief Save background into FITS file
399 *
400 * @param[in] filename Background table FITS file name.
401 * @param[in] clobber Overwrite existing file?
402 *
403 * Saves background into a FITS file. If a file with the given @p filename
404 * does not yet exist it will be created, otherwise the method opens the
405 * existing file. The method will create a (or replace an existing)
406 * background extension. The extension name can be specified as part
407 * of the @p filename, or if no extension name is given, is assumed to be
408 * `BKG`.
409 *
410 * An existing file will only be modified if the @p clobber flag is set to
411 * true.
412 ***************************************************************************/
413void GCTABackground2D::save(const GFilename& filename, const bool& clobber) const
414{
415 // Get extension name
417
418 // Open or create FITS file (without extension name since the requested
419 // extension may not yet exist in the file)
420 GFits fits(filename.url(), true);
421
422 // Remove extension if it exists already
423 if (fits.contains(extname)) {
424 fits.remove(extname);
425 }
426
427 // Create binary table
429
430 // Write the background table
431 write(table);
432
433 // Set binary table extension name
434 table.extname(extname);
435
436 // Append table to FITS file
437 fits.append(table);
438
439 // Save to file
440 fits.save(clobber);
441
442 // Return
443 return;
444}
445
446
447/***********************************************************************//**
448 * @brief Returns MC instrument direction
449 *
450 * @param[in] energy Event energy.
451 * @param[in] time Event trigger time.
452 * @param[in,out] ran Random number generator.
453 * @return CTA instrument direction.
454 *
455 * Returns a Monte Carlo simulated instrument direction for a given @p energy
456 * and event trigger @p time. The simulation is done using a rejection
457 * method that makes use of the background rate access operator. This assures
458 * that the simulation is consistent with the background rate interpolation
459 * that is done by the access operator.
460 ***************************************************************************/
462 const GTime& time,
463 GRan& ran) const
464{
465 // If Monte Carlo has not been initialised then initialise the cache now
466 if (m_mc_max.empty()) {
468 }
469
470 // Determine energy range of response table
471 GEnergy emin = m_mc_spectrum.energy(0);
473 if (energy < emin || energy > emax) {
474 std::string msg = "Event energy "+energy.print()+" is outside the "
475 "energy range ["+emin.print()+", "+emax.print()+"] "
476 "covered by the background response table. Please "
477 "restrict the energy range of the simulation to the "
478 "validity range of the background response table.";
480 }
481
482 // Allocate instrument direction
483 GCTAInstDir dir;
484
485 // Continue only if there is a background model
486 if (!m_mc_max.empty()) {
487
488 // Get reference to node array for the energy axis
489 const GNodeArray& energy_nodes = m_background.axis_nodes(m_inx_energy);
490
491 // Get log10(energy) in TeV
492 double logE = energy.log10TeV();
493
494 // Set values for node arrays
495 energy_nodes.set_value(logE);
496
497 // Get indices for energy interpolation
498 int inx_left = energy_nodes.inx_left();
499 int inx_right = energy_nodes.inx_right();
500
501 // Get maximum background rate as the maximum of the left and the right
502 // energy node. Add some margin.
503 double max_rate = m_mc_max[inx_left];
504 if (m_mc_max[inx_right] > max_rate) {
505 max_rate = m_mc_max[inx_right];
506 }
507 max_rate *= 1.5;
508
509 // Get theta width
510 double theta_width = 2.0 * m_theta_max;
511
512 // Get instrument direction using a rejection method. This assures
513 // that the Monte Carlo sample follows the model distribution
514 while (true) {
515
516 // Get randomized detx and dety in radians
517 double detx = ran.uniform() * theta_width - m_theta_max;
518 double dety = ran.uniform() * theta_width - m_theta_max;
519
520 // Get background rate for these coordinates. If the background
521 // rate is zero then get other coordinates.
522 double value = this->operator()(logE, detx, dety);
523 if (value <= 0.0) {
524 continue;
525 }
526
527 // Dump a warning if value is larger than the maximum rate. This
528 // should never happen!
529 if (value > max_rate) {
530 std::string msg = "Background rate "+gammalib::str(value)+" "
531 "for "+energy.print()+" at "
532 "DETXY=["+
535 "is larger than the maximum expected rate "+
536 gammalib::str(max_rate)+". Something is "
537 "wrong with the code.";
539 }
540
541 // Get uniform random number
542 double uniform = ran.uniform() * max_rate;
543
544 // Exit loop if we're not larger than the background rate value
545 if (uniform <= value) {
546 dir.detx(detx);
547 dir.dety(dety);
548 break;
549 }
550
551 } // endwhile: loop until instrument direction was accepted
552
553 } // endif: there were cube pixels and maps
554
555 // Return instrument direction
556 return dir;
557}
558
559
560/***********************************************************************//**
561 * @brief Returns background rate integrated over energy interval in units
562 * of events s\f$^{-1}\f$ sr\f$^{-1}\f$
563 *
564 * @param[in] dir Instrument direction.
565 * @param[in] emin Minimum energy of energy interval.
566 * @param[in] emax Maximum energy of energy interval.
567 * @return Integrated background count rate (events s\f$^{-1}\f$ sr\f$^{-1}\f$).
568 *
569 * Returns the background count rate for a given instrument direction that
570 * is integrated over a specified energy interval. The rate is given per
571 * ontime.
572 *
573 * If the energy interval is not positive, a zero background rate is
574 * returned.
575 ***************************************************************************/
577 const GEnergy& emin,
578 const GEnergy& emax) const
579{
580 // Initialise rate
581 double rate = 0.0;
582
583 // Continue only if background model is valid and the energy interval is
584 // positive
585 if (is_valid() && (emax > emin)) {
586
587 // Compute theta in radians
588 double theta = std::sqrt(dir.detx() * dir.detx() + dir.dety() * dir.dety());
589
590 // Continue only if theta is within validity range
591 if ((theta >= m_theta_min) && (theta <= m_theta_max)) {
592
593 // Initialise first and second node
594 double x1 = emin.MeV();
595 double x2 = 0.0;
596 double f1 = (*this)(emin.log10TeV(), dir.detx(), dir.dety());
597 double f2 = 0.0;
598
599 // Loop over all nodes
600 for (int i = 0; i < m_energy.size(); ++i) {
601
602 // If node energy is below x1 then skip node
603 if (m_energy[i].MeV() <= x1) {
604 continue;
605 }
606
607 // If node energy is above emax then use emax as energy
608 if (m_energy[i] > emax) {
609 x2 = emax.MeV();
610 f2 = (*this)(emax.log10TeV(), dir.detx(), dir.dety());
611 }
612
613 // ... otherwise use node energy
614 else {
615 x2 = m_energy[i].MeV();
616 f2 = this->rate(i, theta);
617 }
618
619 // Compute integral
620 if (f1 > 0.0 && f2 > 0.0) {
621 rate += gammalib::plaw_integral(x1, f1, x2, f2);
622 }
623
624 // Set second node as first node
625 x1 = x2;
626 f1 = f2;
627
628 // If node energy is above emax then break now
629 if (m_energy[i] > emax) {
630 break;
631 }
632
633 } // endfor: looped over all nodes
634
635 // If last node energy is below emax then compute last part of
636 // integral up to emax
637 if (x1 < emax.MeV()) {
638 x2 = emax.MeV();
639 f2 = (*this)(emax.log10TeV(), dir.detx(), dir.dety());
640 if (f1 > 0.0 && f2 > 0.0) {
641 rate += gammalib::plaw_integral(x1, f1, x2, f2);
642 }
643 }
644
645 } // endif: theta was in valid range
646
647 } // endif: background model was valid and energy interval was positive
648
649 // Return background rate
650 return rate;
651}
652
653
654/***********************************************************************//**
655 * @brief Print background information
656 *
657 * @param[in] chatter Chattiness.
658 * @return String containing background information.
659 ***************************************************************************/
660std::string GCTABackground2D::print(const GChatter& chatter) const
661{
662 // Initialise result string
663 std::string result;
664
665 // Continue only if chatter is not silent
666 if (chatter != SILENT) {
667
668 // Append header
669 result.append("=== GCTABackground2D ===");
670
671 // Append information
672 result.append("\n"+gammalib::parformat("Filename")+m_filename);
673
674 // If background model is valid then print information
675 if (is_valid()) {
676
677 // Compute THETA boundaries in deg
678 double theta_min = m_background.axis_lo(m_inx_theta, 0);
679 double theta_max = m_background.axis_hi(m_inx_theta, m_num_theta-1);
680
681 // Compute energy boundaries in TeV
682 double emin = m_background.axis_lo(m_inx_energy, 0);
684
685 // Append information
686 result.append("\n"+gammalib::parformat("Number of THETA bins") +
688 result.append("\n"+gammalib::parformat("Number of energy bins") +
690 result.append("\n"+gammalib::parformat("Energy range"));
691 result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
692 result.append("\n"+gammalib::parformat("Offset angle range"));
693 result.append(gammalib::str(theta_min)+" - "+gammalib::str(theta_max)+" deg");
694
695 } // endif: there were 2 axis
696
697 // ... otherwise show empty array
698 else {
699 result.append("\n"+gammalib::parformat("Number of THETA bins") +
700 gammalib::str(0));
701 result.append("\n"+gammalib::parformat("Number of energy bins") +
702 gammalib::str(0));
703 }
704
705 } // endif: chatter was not silent
706
707 // Return result
708 return result;
709}
710
711
712/*==========================================================================
713 = =
714 = Private methods =
715 = =
716 ==========================================================================*/
717
718/***********************************************************************//**
719 * @brief Initialise class members
720 ***************************************************************************/
722{
723 // Initialise members
726 m_energy.clear();
727 m_inx_theta = 0;
728 m_inx_energy = 1;
729 m_inx_bgd = 0;
730 m_num_theta = 0;
731 m_num_energy = 0;
732 m_num[0] = 0;
733 m_num[1] = 0;
734 m_theta_min = 0.0;
735 m_theta_max = 0.0;
736 m_logE_min = 0.0;
737 m_logE_max = 0.0;
738
739 // Initialise MC cache
740 m_mc_max.clear();
742
743 // Return
744 return;
745}
746
747
748/***********************************************************************//**
749 * @brief Copy class members
750 *
751 * @param[in] bgd Background.
752 ***************************************************************************/
754{
755 // Copy members
758 m_energy = bgd.m_energy;
761 m_inx_bgd = bgd.m_inx_bgd;
764 m_num[0] = bgd.m_num[0];
765 m_num[1] = bgd.m_num[1];
770
771 // Copy MC cache
772 m_mc_max = bgd.m_mc_max;
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 members from background table
792 *
793 * @exception GException::invalid_value
794 * Response table is not three-dimensional.
795 *
796 * Set class members based on the background table. The THETA axis
797 * will be set to radians, the energy axis will be set to log10.
798 ***************************************************************************/
800{
801 // Throw an exception if the table is not two-dimensional
802 if (m_background.axes() != 2) {
803 std::string msg = "Expected two-dimensional background "
804 "response table but found "+
806 " dimensions. Please specify a two-dimensional "
807 "background.";
809 }
810
811 // Get mandatory indices (throw exception if not found)
812 m_inx_theta = m_background.axis("THETA");
813 m_inx_energy = m_background.axis("ENERG");
815
816 // Get axes dimensions
819
820 // Set dimension array
823
824 // Set THETA axis to radians
826
827 // Set energy axis to logarithmic scale
829
830 // Compute THETA boundaries in radians
835
836 // Compute energy boundaries in log10(TeV)
841 m_logE_min = emin.log10TeV();
842 m_logE_max = emax.log10TeV();
843
844 // Setup energy vector
845 m_energy.clear();
846 for (int i = 0; i < m_num_energy; ++i) {
851 GEnergy energy(std::sqrt(elo.TeV()*ehi.TeV()), "TeV");
852 m_energy.append(energy);
853 }
854
855 // Return
856 return;
857}
858
859
860/***********************************************************************//**
861 * @brief Return background rate bin index
862 *
863 * @param[in] itheta THETA index.
864 * @param[in] iebin Energy index.
865 * @return Background rate bin index.
866 *
867 * Returns the background rate bin index independent of the ordering.
868 ***************************************************************************/
869int GCTABackground2D::index(const int& itheta,
870 const int& iebin) const
871{
872 // Set index arrays
873 int inx[2];
874 inx[m_inx_theta] = itheta;
875 inx[m_inx_energy] = iebin;
876
877 // Compute index
878 int index = inx[0] + inx[1] * m_num[0];
879
880 // Return index
881 return index;
882}
883
884
885/***********************************************************************//**
886 * @brief Initialise Monte Carlo cache
887 *
888 * @exception GException::invalid_value
889 * Background response table was empty.
890 *
891 * Initialises the cache for Monte Carlo sampling.
892 ***************************************************************************/
894{
895 // Set number of subbins per energy bin of response cube
896 const int nsubbins = 10;
897
898 // Initialise cache
899 m_mc_max.clear();
901
902 // Initialise maximum rate array
904
905 // Compute DETX and DETY binsize in radians
906 double detx_bin = 2.0 * m_theta_max / m_num_theta;
907 double dety_bin = 2.0 * m_theta_max / m_num_theta;
908
909 // Set energy bins for Monte Carlo cache spectrum. Make sure that the
910 // actual energy nodes are part of the energy bins to avoid that sharp
911 // structures in the background rate variations are missed.
912 GEnergies energies;
913 for (int i = 0; i < m_num_energy; ++i) {
914
915 // If this is the first node then add energies below the first node
916 // down to the lower energy limit of the background response
917 if (i == 0) {
920 GEnergies enodes(nsubbins/2, emin, m_energy[i]);
921 enodes.remove(enodes.size()-1); // avoids duplication of energy
922 energies.extend(enodes);
923 }
924
925 // If this is not the last node then append nodes between current
926 // node and next node
927 if (i < m_num_energy-1) {
928 GEnergies enodes(nsubbins, m_energy[i], m_energy[i+1]);
929 enodes.remove(enodes.size()-1); // avoids duplication of energy
930 energies.extend(enodes);
931 }
932
933 // ... otherwise this is the last node. Add energies above the last
934 // node up to the upper energy limit of the background response
935 else {
938 GEnergies enodes(nsubbins/2, m_energy[i], emax);
939 energies.extend(enodes);
940 }
941
942 } // endfor: looped over energy bins of background response
943
944 // Loop over all energy bins
945 for (int i = 0; i < energies.size(); ++i) {
946
947 // Get logE of energy bin
948 double logE = energies[i].log10TeV();
949
950 // Initialise cache with cumulative pixel fluxes and compute total flux
951 // in response table for normalization. Negative pixels are excluded
952 // from the cumulative map.
953 double total_flux = 0.0; // units: events/s/MeV
954 double dx = 0.5 * detx_bin;
955 double dy = 0.5 * dety_bin;
956 double detx = -m_theta_max + dx;
957 for (int ix = 0; ix < 2*m_num_theta; ++ix, detx += detx_bin) {
958 double dety = -m_theta_max + dy;
959 for (int iy = 0; iy < 2*m_num_theta; ++iy, dety += dety_bin) {
960
961 // Compute intensities
962 double i0 = (*this)(logE, detx, dety);
963 double i1 = (*this)(logE, detx-dx, dety-dy);
964 double i2 = (*this)(logE, detx, dety-dy);
965 double i3 = (*this)(logE, detx+dx, dety-dy);
966 double i4 = (*this)(logE, detx+dx, dety);
967 double i5 = (*this)(logE, detx+dx, dety+dy);
968 double i6 = (*this)(logE, detx, dety+dy);
969 double i7 = (*this)(logE, detx-dx, dety+dy);
970 double i8 = (*this)(logE, detx-dx, dety);
971
972 // Compute solid angle of the 8 pixel wedges
973 double s1 = solid_angle(detx, dety, detx-dx, dety-dy, detx, dety-dy);
974 double s2 = solid_angle(detx, dety, detx, dety-dy, detx+dx, dety-dy);
975 double s3 = solid_angle(detx, dety, detx+dx, dety-dy, detx+dx, dety);
976 double s4 = solid_angle(detx, dety, detx+dx, dety, detx+dx, dety+dy);
977 double s5 = solid_angle(detx, dety, detx+dx, dety+dy, detx, dety+dy);
978 double s6 = solid_angle(detx, dety, detx, dety+dy, detx-dx, dety+dy);
979 double s7 = solid_angle(detx, dety, detx-dx, dety+dy, detx-dx, dety);
980 double s8 = solid_angle(detx, dety, detx-dx, dety, detx-dx, dety-dy);
981
982 // Compute flux by summing the flux in 8 pixel wedges
983 double flux1 = s1 * (i1 + i2 + i0);
984 double flux2 = s2 * (i2 + i3 + i0);
985 double flux3 = s3 * (i3 + i4 + i0);
986 double flux4 = s4 * (i4 + i5 + i0);
987 double flux5 = s5 * (i5 + i6 + i0);
988 double flux6 = s6 * (i6 + i7 + i0);
989 double flux7 = s7 * (i7 + i8 + i0);
990 double flux8 = s8 * (i8 + i1 + i0);
991 double flux = (flux1 + flux2 + flux3 + flux4 +
992 flux5 + flux6 + flux7 + flux8) / 3.0;
993
994 // Sum flux
995 if (flux > 0.0) {
996 total_flux += flux;
997 }
998
999 } // endfor: loop over DETY
1000
1001 } // endfor: loop over DETX
1002
1003 // Append node if flux is positive
1004 if (total_flux > 0.0) {
1005 m_mc_spectrum.append(energies[i], total_flux);
1006 }
1007
1008 // Dump spectrum for debugging
1009 #if defined(G_DEBUG_MC_INIT)
1010 std::cout << "Energy=" << energies[i];
1011 std::cout << " Rate_total=" << total_flux;
1012 std::cout << " events/(s MeV)" << std::endl;
1013 #endif
1014
1015 } // endfor: looped over energy bins
1016
1017 // Make sure that spectrum is not empty (if this is the case the entire
1018 // background response table was empty)
1019 if (m_mc_spectrum.nodes() == 0) {
1020 std::string msg = "Background response table is empty. Please provide "
1021 "a valid three-dimensional background template.";
1023 }
1024
1025 // Return
1026 return;
1027}
1028
1029
1030/***********************************************************************//**
1031 * @brief Initialise array of maximum background rate
1032 *
1033 * Initialises the array m_mc_max that holds the maximum background rate for
1034 * each energy of the background model.
1035 ***************************************************************************/
1037{
1038 // Initialise cache
1039 m_mc_max.clear();
1040
1041 // Loop over all maps to determine the maximum map intensity in units
1042 // of events/s/sr/MeV
1043 for (int iebin = 0; iebin < m_num_energy; ++iebin) {
1044
1045 // Initialise maximum rate
1046 double max_rate = 0.0;
1047
1048 // Loop over THETA angles
1049 for (int itheta = 0; itheta < m_num_theta; ++itheta) {
1050
1051 // Get bin index
1052 int inx = index(itheta, iebin);
1053
1054 // Get background rate
1055 double rate = m_background(m_inx_bgd, inx);
1056
1057 // If background rate is larger than maximum then store the
1058 // background rate as new maximum
1059 if (rate > max_rate) {
1060 max_rate = rate;
1061 }
1062
1063 } // endfor: looped over THETA
1064
1065 // Append maximum rate
1066 m_mc_max.push_back(max_rate);
1067
1068 } // endif: looped over all maps
1069
1070 // Return
1071 return;
1072}
1073
1074
1075/***********************************************************************//**
1076 * @brief Compute solid angle of pixel wedge
1077 *
1078 * @param[in] detx1 DETX of first edge.
1079 * @param[in] dety1 DETY of first edge.
1080 * @param[in] detx2 DETX of second edge.
1081 * @param[in] dety2 DETY of second edge.
1082 * @param[in] detx3 DETX of third edge.
1083 * @param[in] dety3 DETY of third edge.
1084 * @return Solid angle (steradians).
1085 *
1086 * Estimate the solid angle subtended by 3 coordinates using Huilier's
1087 * theorem.
1088 *
1089 * Below, the definiton of the pixel cornes and sides are shown as used
1090 * within the code.
1091 *
1092 * a12
1093 * 1---------2
1094 * | /
1095 * | /
1096 * | /
1097 * | /
1098 * a13| /a23
1099 * | /
1100 * | /
1101 * | /
1102 * |/
1103 * 3
1104 *
1105 ***************************************************************************/
1106double GCTABackground2D::solid_angle(const double& detx1, const double& dety1,
1107 const double& detx2, const double& dety2,
1108 const double& detx3, const double& dety3) const
1109{
1110 // Set sky directions
1111 double theta1 = std::sqrt(detx1 * detx1 + dety1 * dety1);
1112 double phi1 = std::atan2(dety1, detx1);
1113 double theta2 = std::sqrt(detx2 * detx2 + dety2 * dety2);
1114 double phi2 = std::atan2(dety2, detx2);
1115 double theta3 = std::sqrt(detx3 * detx3 + dety3 * dety3);
1116 double phi3 = std::atan2(dety3, detx3);
1117 GSkyDir dir1;
1118 GSkyDir dir2;
1119 GSkyDir dir3;
1120 dir1.radec(phi1, gammalib::pihalf - theta1);
1121 dir2.radec(phi2, gammalib::pihalf - theta2);
1122 dir3.radec(phi3, gammalib::pihalf - theta3);
1123
1124 // Compute angular distances between pixel corners
1125 double a12 = dir1.dist(dir2);
1126 double a13 = dir1.dist(dir3);
1127 double a23 = dir2.dist(dir3);
1128
1129 // Compute solid angle
1130 double s = 0.5 * (a12 + a23 + a13);
1131 double solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
1132 std::tan(0.5*(s-a12)) *
1133 std::tan(0.5*(s-a23)) *
1134 std::tan(0.5*(s-a13))));
1135
1136 // Return solid angle
1137 return solidangle;
1138}
1139
1140
1141/***********************************************************************//**
1142 * @brief Return background rate for a given energy bin and offset angle
1143 * value (events/s/MeV/sr)
1144 *
1145 * @param[in] iebin Energy index.
1146 * @param[in] theta Offset angle in nominal system (radians).
1147 * @return Background rate (events/s/MeV/sr)
1148 ***************************************************************************/
1149double GCTABackground2D::rate(const int& iebin, const double& theta) const
1150{
1151 // Retrieve reference to node array
1153
1154 // Set offset angle value for node array interpolation
1155 thetas.set_value(theta);
1156
1157 // Get bin indices
1158 int inx_left = index(thetas.inx_left(), iebin);
1159 int inx_right = index(thetas.inx_right(), iebin);
1160
1161 // Bi-linearly interpolate the rates
1162 double rate = thetas.wgt_left() * m_background(m_inx_bgd, inx_left) +
1163 thetas.wgt_right() * m_background(m_inx_bgd, inx_right);
1164
1165 // Make sure that background rate is not negative
1166 if (rate < 0.0) {
1167 rate = 0.0;
1168 }
1169
1170 // Return
1171 return rate;
1172}
#define G_INIT_MC_CACHE
CTA 2D background class definition.
CTA instrument direction class interface definition.
Definition of support function used by CTA classes.
Exception handler interface definition.
#define G_SET_MEMBERS
Filename class interface definition.
FITS binary table class definition.
FITS file class interface definition.
Mathematical function definitions.
Random number generator class definition.
Gammalib tools definition.
GChatter
Definition GTypemaps.hpp:33
@ SILENT
Definition GTypemaps.hpp:34
CTA 2D background class.
double rate(const int &iebin, const double &theta) const
Return background rate for a given energy bin and offset angle value (events/s/MeV/sr)
void save(const GFilename &filename, const bool &clobber=false) const
Save background into FITS file.
std::string print(const GChatter &chatter=NORMAL) const
Print background information.
const GCTAResponseTable & table(void) const
Return response table.
int index(const int &itheta, const int &iebin) const
Return background rate bin index.
virtual ~GCTABackground2D(void)
Destructor.
GFilename filename(void) const
Return filename.
int m_inx_theta
THETA axis index.
bool is_valid(void) const
Return validity of background model.
GCTABackground2D(void)
Void constructor.
GFilename m_filename
Name of background response file.
void init_mc_max_rate(void) const
Initialise array of maximum background rate.
void init_members(void)
Initialise class members.
void clear(void)
Clear background.
GCTAResponseTable m_background
Background response table.
void read(const GFitsTable &table)
Read background from FITS table.
double m_theta_min
THETA minimum (radians)
void free_members(void)
Delete class members.
void write(GFitsBinTable &table) const
Write background into FITS table.
void copy_members(const GCTABackground2D &bgd)
Copy class members.
GEnergies m_energy
Vector of energies.
int m_num_energy
Number of energy bins.
double m_logE_min
Log10(E/TeV) minimum.
void init_mc_cache(void) const
Initialise Monte Carlo cache.
double m_theta_max
THETA maximum (radians)
GCTABackground2D & operator=(const GCTABackground2D &bgd)
Assignment operator.
GCTAInstDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Returns MC instrument direction.
int m_num[2]
Array of number of bins.
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.
GCTABackground2D * clone(void) const
Clone background.
int m_inx_bgd
Background index.
void set_members(void)
Set members from background table.
int m_num_theta
Number of THETA bins.
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 .
int m_inx_energy
Energy axis index.
double m_logE_max
Log10(E/TeV) maximum.
std::vector< double > m_mc_max
Maximum background rate.
GModelSpectralNodes m_mc_spectrum
Response cube spectrum.
virtual double operator()(const double &logE, const double &detx, const double &dety) const
Return background rate in units of events MeV s sr .
void load(const GFilename &filename)
Load background from FITS file.
Abstract base class for the CTA background model.
void init_members(void)
Initialise class members.
GCTABackground & operator=(const GCTABackground &bgd)
Assignment operator.
void free_members(void)
Delete class members.
CTA instrument direction class.
void dety(const double &y)
Set DETY coordinate (in radians)
void detx(const double &x)
Set DETX coordinate (in radians)
GCTAResponseTable * clone(void) const
Clone response table.
void read(const GFitsTable &table)
Read response table from FITS table HDU.
const std::string & axis_lo_unit(const int &axis) const
Return lower bin boundary unit for axis.
int table(const std::string &name) const
Determine index of table.
const std::string & axis_hi_unit(const int &axis) const
Return upper bin boundary unit for axis.
void axis_log10(const int &axis)
Set nodes for a logarithmic (base 10) axis.
void axis_radians(const int &axis)
Set nodes for a radians axis.
int axis(const std::string &name) const
Determine index of an axis.
const int & axes(void) const
Return number of axes of the tables.
void write(GFitsTable &table) const
Write response table into FITS table HDU.
int axis_bins(const int &axis) const
Return number bins in an axis.
const double & axis_hi(const int &axis, const int &bin) const
Return upper bin boundary for bin in axis.
const GNodeArray & axis_nodes(const int &axis) const
Return axis nodes.
void clear(void)
Clear response table.
const double & axis_lo(const int &axis, const int &bin) const
Return lower bin boundary for bin in axis.
Energy container class.
Definition GEnergies.hpp:60
void clear(void)
Clear energy container.
int size(void) const
Return number of energies in container.
void extend(const GEnergies &energies)
Append energy container.
GEnergy & append(const GEnergy &energy)
Append energy to container.
void remove(const int &index)
Remove energy from container.
Class that handles energies in a unit independent way.
Definition GEnergy.hpp:48
double log10TeV(void) const
Return log10 of energy in TeV.
Definition GEnergy.cpp:457
double MeV(void) const
Return energy in MeV.
Definition GEnergy.cpp:321
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition GEnergy.cpp:748
double TeV(void) const
Return energy in TeV.
Definition GEnergy.cpp:348
Filename class.
Definition GFilename.hpp:62
std::string url(void) const
Return Uniform Resource Locator (URL)
std::string extname(const std::string &defaultname="") const
Return extension name.
void clear(void)
Clear file name.
FITS binary table class.
Abstract interface for FITS table.
void remove(const int &colnum)
Remove column from the table.
FITS file class.
Definition GFits.hpp:63
bool contains(const int &extno) const
Check if HDU exists in FITS file.
Definition GFits.hpp:282
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition GFits.cpp:678
void close(void)
Close FITS file.
Definition GFits.cpp:1342
void remove(const int &extno)
Remove HDU from FITS file.
Definition GFits.cpp:862
void save(const bool &clobber=false)
Saves FITS file.
Definition GFits.cpp:1178
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition GFits.cpp:482
virtual void clear(void)
Clear spectral nodes model.
int nodes(void) const
Return number of nodes.
void append(const GEnergy &energy, const double &intensity)
Append node.
GEnergy energy(const int &index) const
Return node energy.
Node array class.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
const int & inx_right(void) const
Returns right node index.
const int & inx_left(void) const
Returns left node index.
const double & wgt_right(void) const
Returns right node weight.
const double & wgt_left(void) const
Returns left node weight.
Random number generator class.
Definition GRan.hpp:44
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition GRan.cpp:242
Sky direction class.
Definition GSkyDir.hpp:62
void radec(const double &ra, const double &dec)
Set equatorial sky direction (radians)
Definition GSkyDir.cpp:197
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition GSkyDir.hpp:271
Time class.
Definition GTime.hpp:55
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition GTools.cpp:1143
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 str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition GTools.cpp:489
const double pihalf
Definition GMath.hpp:38
const double deg2rad
Definition GMath.hpp:43
const std::string extname_cta_background2d
std::string gadf_hduclas4(const GFits &fits, const std::string &hduclas4)
Return extension name for GADF response table of given HDU class 4.
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition GTools.cpp:1386
const double rad2deg
Definition GMath.hpp:44