Justine Devin / gammalib

Commit 84df48dd5242131ecf9637d669f7df605f4d80ee

Authored by Jürgen Knödlseder
1 parent 1f5e1957

Implement background model support in GCTAResponse and GCTAObservation. 

The CTA background model, represented by the abstract base class GCTABackground, has been added as 4th component of the instrument response. The CTA observation definition file has been updated by adding a "Background" parameter to the observation tag, allowing the specification of a filename for the background model.
Inline Side-by-side
Showing 9 changed files with 174 additions and 75 deletions
doc/source/user_manual/modules/cta.rst
  View file @84df48d
... ... @@ -98,6 +98,7 @@ the format:
98 98 <parameter name="EffectiveArea" file="..."/>
99 99 <parameter name="PointSpreadFunction" file="..."/>
100 100 <parameter name="EnergyDispersion" file="..."/>
  101 + <parameter name="Background" file="..."/>
101 102 </observation>
102 103  
103 104 for an unbinned observation and:
... ... @@ -109,6 +110,7 @@ for an unbinned observation and:
109 110 <parameter name="EffectiveArea" file="..."/>
110 111 <parameter name="PointSpreadFunction" file="..."/>
111 112 <parameter name="EnergyDispersion" file="..."/>
  113 + <parameter name="Background" file="..."/>
112 114 </observation>
113 115  
114 116 for a binned observation. Each parameter within the ``<observation>`` tag
... ... @@ -129,6 +131,8 @@ about the off-axis reponse and event selection (see
129 131 So far, energy dispersion is not supported by GammaLib, hence the
130 132 ``file`` attribute for the ``EnergyDispersion`` parameter can be left
131 133 blank.
  134 +The ``Background`` parameter allows specifying of a file that contains
  135 +the instrument background model for the observation.
132 136  
133 137 The ``instrument`` attribute of CTA observations can be one of
134 138 ``CTA``, ``HESS``, ``MAGIC`` or ``VERITAS``. This allows mixing of
... ...
inst/cta/include/GCTAResponse.hpp
  View file @84df48d
... ... @@ -51,6 +51,7 @@ class GCTAInstDir;
51 51 class GCTAAeff;
52 52 class GCTAPsf;
53 53 class GCTAEdisp;
  54 +class GCTABackground;
54 55  
55 56  
56 57 /***********************************************************************//**
... ... @@ -103,24 +104,27 @@ public:
103 104 virtual GEbounds ebounds_src(const GEnergy& obsEnergy) const;
104 105  
105 106 // Other Methods
106   - GCTAEventAtom* mc(const double& area, const GPhoton& photon,
107   - const GObservation& obs, GRan& ran) const;
108   - void caldb(const GCaldb& caldb);
109   - const GCaldb& caldb(void) const;
110   - void load(const std::string& rspname);
111   - void eps(const double& eps);
112   - const double& eps(void) const;
113   - void load_aeff(const std::string& filename);
114   - void load_psf(const std::string& filename);
115   - void load_edisp(const std::string& filename);
116   - void offset_sigma(const double& sigma);
117   - double offset_sigma(void) const;
118   - const GCTAAeff* aeff(void) const;
119   - void aeff(GCTAAeff* aeff);
120   - const GCTAPsf* psf(void) const;
121   - void psf(GCTAPsf* psf);
122   - const GCTAEdisp* edisp(void) const;
123   - void edisp(GCTAEdisp* edisp);
  107 + GCTAEventAtom* mc(const double& area, const GPhoton& photon,
  108 + const GObservation& obs, GRan& ran) const;
  109 + void caldb(const GCaldb& caldb);
  110 + const GCaldb& caldb(void) const;
  111 + void load(const std::string& rspname);
  112 + void eps(const double& eps);
  113 + const double& eps(void) const;
  114 + void load_aeff(const std::string& filename);
  115 + void load_psf(const std::string& filename);
  116 + void load_edisp(const std::string& filename);
  117 + void load_background(const std::string& filename);
  118 + void offset_sigma(const double& sigma);
  119 + double offset_sigma(void) const;
  120 + const GCTAAeff* aeff(void) const;
  121 + void aeff(GCTAAeff* aeff);
  122 + const GCTAPsf* psf(void) const;
  123 + void psf(GCTAPsf* psf);
  124 + const GCTAEdisp* edisp(void) const;
  125 + void edisp(GCTAEdisp* edisp);
  126 + const GCTABackground* background(void) const;
  127 + void background(GCTABackground* background);
124 128  
125 129 // Low-level response methods
126 130 double aeff(const double& theta,
... ... @@ -172,13 +176,14 @@ private:
172 176 std::string irf_filename(const std::string& filename) const;
173 177  
174 178 // Private data members
175   - GCaldb m_caldb; //!< Calibration database
176   - std::string m_rspname; //!< Name of the instrument response
177   - double m_eps; //!< Integration precision
178   - GCTAAeff* m_aeff; //!< Effective area
179   - GCTAPsf* m_psf; //!< Point spread function
180   - GCTAEdisp* m_edisp; //!< Energy dispersion
181   - mutable bool m_apply_edisp; //!< Apply energy dispersion
  179 + GCaldb m_caldb; //!< Calibration database
  180 + std::string m_rspname; //!< Name of the instrument response
  181 + double m_eps; //!< Integration precision
  182 + GCTAAeff* m_aeff; //!< Effective area
  183 + GCTAPsf* m_psf; //!< Point spread function
  184 + GCTAEdisp* m_edisp; //!< Energy dispersion
  185 + GCTABackground* m_background; //!< Energy dispersion
  186 + mutable bool m_apply_edisp; //!< Apply energy dispersion
182 187  
183 188 // Npred cache
184 189 mutable std::vector<std::string> m_npred_names; //!< Model names
... ... @@ -366,4 +371,31 @@ void GCTAResponse::edisp(GCTAEdisp* edisp)
366 371 return;
367 372 }
368 373  
  374 +
  375 +
  376 +
  377 +/***********************************************************************//**
  378 + * @brief Return pointer to background model
  379 + *
  380 + * @return Pointer to background model.
  381 + ***************************************************************************/
  382 +inline
  383 +const GCTABackground* GCTAResponse::background(void) const
  384 +{
  385 + return m_background;
  386 +}
  387 +
  388 +
  389 +/***********************************************************************//**
  390 + * @brief Set pointer to background model
  391 + *
  392 + * @param[in] edisp Pointer to background model.
  393 + ***************************************************************************/
  394 +inline
  395 +void GCTAResponse::background(GCTABackground* background)
  396 +{
  397 + m_background = background;
  398 + return;
  399 +}
  400 +
369 401 #endif /* GCTARESPONSE_HPP */
... ...
inst/cta/pyext/GCTAResponse.i
  View file @84df48d
... ... @@ -74,24 +74,27 @@ public:
74 74 virtual GEbounds ebounds_src(const GEnergy& obsEnergy) const;
75 75  
76 76 // Other Methods
77   - GCTAEventAtom* mc(const double& area, const GPhoton& photon,
78   - const GObservation& obs, GRan& ran) const;
79   - void caldb(const GCaldb& caldb);
80   - const GCaldb& caldb(void) const;
81   - void load(const std::string& rspname);
82   - void eps(const double& eps);
83   - const double& eps(void) const;
84   - void load_aeff(const std::string& filename);
85   - void load_psf(const std::string& filename);
86   - void load_edisp(const std::string& filename);
87   - void offset_sigma(const double& sigma);
88   - double offset_sigma(void) const;
89   - const GCTAAeff* aeff(void) const;
90   - void aeff(GCTAAeff* aeff);
91   - const GCTAPsf* psf(void) const;
92   - void psf(GCTAPsf* psf);
93   - const GCTAEdisp* edisp(void) const;
94   - void edisp(GCTAEdisp* edisp);
  77 + GCTAEventAtom* mc(const double& area, const GPhoton& photon,
  78 + const GObservation& obs, GRan& ran) const;
  79 + void caldb(const GCaldb& caldb);
  80 + const GCaldb& caldb(void) const;
  81 + void load(const std::string& rspname);
  82 + void eps(const double& eps);
  83 + const double& eps(void) const;
  84 + void load_aeff(const std::string& filename);
  85 + void load_psf(const std::string& filename);
  86 + void load_edisp(const std::string& filename);
  87 + void load_background(const std::string& filename);
  88 + void offset_sigma(const double& sigma);
  89 + double offset_sigma(void) const;
  90 + const GCTAAeff* aeff(void) const;
  91 + void aeff(GCTAAeff* aeff);
  92 + const GCTAPsf* psf(void) const;
  93 + void psf(GCTAPsf* psf);
  94 + const GCTAEdisp* edisp(void) const;
  95 + void edisp(GCTAEdisp* edisp);
  96 + const GCTABackground* background(void) const;
  97 + void background(GCTABackground* background);
95 98  
96 99 // Low-level response methods
97 100 double aeff(const double& theta,
... ...
inst/cta/src/GCTAObservation.cpp
  View file @84df48d
... ... @@ -44,6 +44,7 @@
44 44 #include "GCTAAeffArf.hpp"
45 45 #include "GCTAAeffPerfTable.hpp"
46 46 #include "GCTAEdisp.hpp"
  47 +#include "GCTABackground.hpp"
47 48  
48 49 /* __ Globals ____________________________________________________________ */
49 50 const GCTAObservation g_obs_cta_seed("CTA");
... ... @@ -288,6 +289,7 @@ void GCTAObservation::response(const std::string&amp; rspname, const GCaldb&amp; caldb)
288 289 * <parameter name="EffectiveArea" file="..."/>
289 290 * <parameter name="PointSpreadFunction" file="..."/>
290 291 * <parameter name="EnergyDispersion" file="..."/>
  292 + * <parameter name="Background" file="..."/>
291 293 * </observation>
292 294 *
293 295 * for an unbinned observation and
... ... @@ -297,12 +299,12 @@ void GCTAObservation::response(const std::string&amp; rspname, const GCaldb&amp; caldb)
297 299 * <parameter name="EffectiveArea" file="..."/>
298 300 * <parameter name="PointSpreadFunction" file="..."/>
299 301 * <parameter name="EnergyDispersion" file="..."/>
  302 + * <parameter name="Background" file="..."/>
300 303 * </observation>
301 304 *
302 305 * for a binned observation.
303 306 *
304 307 * @todo Still supports old ARF, PSF and RMF parameter names.
305   - * @todo PSF correction for ARF not yet implemented.
306 308 ***************************************************************************/
307 309 void GCTAObservation::read(const GXmlElement& xml)
308 310 {
... ... @@ -315,14 +317,14 @@ void GCTAObservation::read(const GXmlElement&amp; xml)
315 317 // Determine number of parameter nodes in XML element
316 318 int npars = xml.elements("parameter");
317 319  
318   - // Verify that XML element has exactly 4 parameters
319   - if (xml.elements() != 4 || npars != 4) {
  320 + // Verify that XML element has exactly 5 parameters
  321 + if (xml.elements() != 5 || npars != 5) {
320 322 throw GException::xml_invalid_parnum(G_READ, xml,
321   - "CTA observation requires exactly 4 parameters.");
  323 + "CTA observation requires exactly 5 parameters.");
322 324 }
323 325  
324 326 // Extract parameters
325   - int npar[] = {0, 0, 0, 0};
  327 + int npar[] = {0, 0, 0, 0, 0};
326 328 for (int i = 0; i < npars; ++i) {
327 329  
328 330 // Get parameter element
... ... @@ -443,13 +445,29 @@ void GCTAObservation::read(const GXmlElement&amp; xml)
443 445 npar[3]++;
444 446 }
445 447  
  448 + // Handle background model
  449 + else if (par->attribute("name") == "Background") {
  450 +
  451 + // Get filename
  452 + std::string filename = par->attribute("file");
  453 +
  454 + // If filename is not empty then load background model
  455 + if (gammalib::strip_whitespace(filename).length() > 0) {
  456 + m_response.load_background(filename);
  457 + }
  458 +
  459 + // Increase number of parameters
  460 + npar[4]++;
  461 + }
  462 +
446 463 } // endfor: looped over all parameters
447 464  
448 465 // Verify that all parameters were found
449   - if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1) {
  466 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1 || npar[4] != 1) {
450 467 throw GException::xml_invalid_parnames(G_READ, xml,
451 468 "Require \"EventList\" or \"CountsMap\" and \"EffectiveArea\""
452   - ", \"PointSpreadFunction\" and \"EnergyDispersion\" parameters.");
  469 + ", \"PointSpreadFunction\", \"EnergyDispersion\" and"
  470 + " \"Background\" parameters.");
453 471 }
454 472  
455 473 // If we have an ARF then remove thetacut if necessary
... ... @@ -485,6 +503,7 @@ void GCTAObservation::read(const GXmlElement&amp; xml)
485 503 * <parameter name="EffectiveArea" file="..."/>
486 504 * <parameter name="PointSpreadFunction" file="..."/>
487 505 * <parameter name="EnergyDispersion" file="..."/>
  506 + * <parameter name="Background" file="..."/>
488 507 * </observation>
489 508 *
490 509 * for an unbinned observation and
... ... @@ -494,6 +513,7 @@ void GCTAObservation::read(const GXmlElement&amp; xml)
494 513 * <parameter name="EffectiveArea" file="..."/>
495 514 * <parameter name="PointSpreadFunction" file="..."/>
496 515 * <parameter name="EnergyDispersion" file="..."/>
  516 + * <parameter name="Background" file="..."/>
497 517 * </observation>
498 518 *
499 519 * for a binned observation.
... ... @@ -514,7 +534,7 @@ void GCTAObservation::write(GXmlElement&amp; xml) const
514 534 // Set event list flag
515 535 bool is_list = (list != NULL);
516 536  
517   - // If XML element has 0 nodes then append 4 parameter nodes
  537 + // If XML element has 0 nodes then append 5 parameter nodes
518 538 if (xml.elements() == 0) {
519 539 if (is_list) {
520 540 xml.append(GXmlElement("parameter name=\"EventList\""));
... ... @@ -525,17 +545,18 @@ void GCTAObservation::write(GXmlElement&amp; xml) const
525 545 xml.append(GXmlElement("parameter name=\"EffectiveArea\""));
526 546 xml.append(GXmlElement("parameter name=\"PointSpreadFunction\""));
527 547 xml.append(GXmlElement("parameter name=\"EnergyDispersion\""));
  548 + xml.append(GXmlElement("parameter name=\"Background\""));
528 549 }
529 550  
530   - // Verify that XML element has exactly 4 parameters
531   - if (xml.elements() != 4 || xml.elements("parameter") != 4) {
  551 + // Verify that XML element has exactly 5 parameters
  552 + if (xml.elements() != 5 || xml.elements("parameter") != 5) {
532 553 throw GException::xml_invalid_parnum(G_WRITE, xml,
533   - "CTA observation requires exactly 4 parameters.");
  554 + "CTA observation requires exactly 5 parameters.");
534 555 }
535 556  
536 557 // Set or update parameter attributes
537   - int npar[] = {0, 0, 0, 0};
538   - for (int i = 0; i < 4; ++i) {
  558 + int npar[] = {0, 0, 0, 0, 0};
  559 + for (int i = 0; i < 5; ++i) {
539 560  
540 561 // Get parameter element
541 562 GXmlElement* par = xml.element("parameter", i);
... ... @@ -619,13 +640,24 @@ void GCTAObservation::write(GXmlElement&amp; xml) const
619 640 npar[3]++;
620 641 }
621 642  
  643 + // Handle Background
  644 + else if (par->attribute("name") == "Background") {
  645 + std::string filename = "";
  646 + if (m_response.background() != NULL) {
  647 + filename = m_response.background()->filename();
  648 + }
  649 + par->attribute("file", filename);
  650 + npar[4]++;
  651 + }
  652 +
622 653 } // endfor: looped over all parameters
623 654  
624 655 // Verify that all required parameters are present
625   - if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1) {
  656 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1 || npar[4] != 1) {
626 657 throw GException::xml_invalid_parnames(G_WRITE, xml,
627   - "Require \"EventList\" or \"CountsMap\" and \"ARF\", \"RMF\""
628   - " \"PSF\" parameters.");
  658 + "Require \"EventList\" or \"CountsMap\" and"
  659 + " \"PointSpreadFunction\", \"EnergyDispersion\""
  660 + " and \"Background\" parameters.");
629 661 }
630 662  
631 663 // Return
... ...
inst/cta/src/GCTAResponse.cpp
  View file @84df48d
... ... @@ -58,6 +58,7 @@
58 58 #include "GCTAPsf.hpp"
59 59 #include "GCTAEdisp.hpp"
60 60 #include "GCTAEdispPerfTable.hpp"
  61 +#include "GCTABackground.hpp"
61 62  
62 63 /* __ Method name definitions ____________________________________________ */
63 64 #define G_CALDB "GCTAResponse::caldb(GCaldb&)"
... ... @@ -737,7 +738,7 @@ void GCTAResponse::load_psf(const std::string&amp; filename)
737 738 ***************************************************************************/
738 739 void GCTAResponse::load_edisp(const std::string& filename)
739 740 {
740   - // Free any existing point spread function instance
  741 + // Free any existing energy dispersion instance
741 742 if (m_edisp != NULL) delete m_edisp;
742 743 m_edisp = NULL;
743 744  
... ... @@ -752,6 +753,25 @@ void GCTAResponse::load_edisp(const std::string&amp; filename)
752 753 m_edisp = new GCTAEdispPerfTable(filename);
753 754 }
754 755  
  756 + // Return
  757 + return;
  758 +}
  759 +
  760 +
  761 +/***********************************************************************//**
  762 + * @brief Load background model
  763 + *
  764 + * @param[in] filename Background model file name.
  765 + *
  766 + * @todo Method to be implemented.
  767 + ***************************************************************************/
  768 +void GCTAResponse::load_background(const std::string& filename)
  769 +{
  770 + // Free any existing background model instance
  771 + if (m_background != NULL) delete m_background;
  772 + m_background = NULL;
  773 +
  774 + // Return
755 775 return;
756 776 }
757 777  
... ... @@ -2349,6 +2369,7 @@ void GCTAResponse::init_members(void)
2349 2369 m_aeff = NULL;
2350 2370 m_psf = NULL;
2351 2371 m_edisp = NULL;
  2372 + m_background = NULL;
2352 2373 m_apply_edisp = false; //!< Switched off by default
2353 2374  
2354 2375 // Initialise Npred cache
... ... @@ -2382,9 +2403,10 @@ void GCTAResponse::copy_members(const GCTAResponse&amp; rsp)
2382 2403 m_npred_values = rsp.m_npred_values;
2383 2404  
2384 2405 // Clone members
2385   - m_aeff = (rsp.m_aeff != NULL) ? rsp.m_aeff->clone() : NULL;
2386   - m_psf = (rsp.m_psf != NULL) ? rsp.m_psf->clone() : NULL;
2387   - m_edisp = (rsp.m_edisp != NULL) ? rsp.m_edisp->clone() : NULL;
  2406 + m_aeff = (rsp.m_aeff != NULL) ? rsp.m_aeff->clone() : NULL;
  2407 + m_psf = (rsp.m_psf != NULL) ? rsp.m_psf->clone() : NULL;
  2408 + m_edisp = (rsp.m_edisp != NULL) ? rsp.m_edisp->clone() : NULL;
  2409 + m_background = (rsp.m_background != NULL) ? rsp.m_background->clone() : NULL;
2388 2410  
2389 2411 // Return
2390 2412 return;
... ... @@ -2397,14 +2419,16 @@ void GCTAResponse::copy_members(const GCTAResponse&amp; rsp)
2397 2419 void GCTAResponse::free_members(void)
2398 2420 {
2399 2421 // Free memory
2400   - if (m_aeff != NULL) delete m_aeff;
2401   - if (m_psf != NULL) delete m_psf;
2402   - if (m_edisp != NULL) delete m_edisp;
  2422 + if (m_aeff != NULL) delete m_aeff;
  2423 + if (m_psf != NULL) delete m_psf;
  2424 + if (m_edisp != NULL) delete m_edisp;
  2425 + if (m_background != NULL) delete m_background;
2403 2426  
2404 2427 // Initialise pointers
2405   - m_aeff = NULL;
2406   - m_psf = NULL;
2407   - m_edisp = NULL;
  2428 + m_aeff = NULL;
  2429 + m_psf = NULL;
  2430 + m_edisp = NULL;
  2431 + m_background = NULL;
2408 2432  
2409 2433 // Return
2410 2434 return;
... ...
inst/cta/test/data/irf_binned.xml
  View file @84df48d
1 1 <observation_list title="observation library">
2 2 <observation name="Crab" id="000001" instrument="CTA">
3 3 <parameter name="CountsMap" file="$(PACKAGE_SOURCE)/inst/cta/test/data/crab_cntmap.fits"/>
4   - <parameter name="EffectiveArea" file="$(PACKAGE_SOURCE)/inst/cta/caldb/dummy/bcf/hard/000001/ea_test.fits"/>
5   - <parameter name="PointSpreadFunction" file="$(PACKAGE_SOURCE)/inst/cta/caldb/dummy/bcf/hard/000001/psf_test.fits"/>
6   - <parameter name="EnergyDispersion" file="$(PACKAGE_SOURCE)/inst/cta/caldb/dummy/bcf/hard/000001/edisp_test.fits"/>
  4 + <parameter name="EffectiveArea" file="$(PACKAGE_SOURCE)/inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h/irf_file.fits"/>
  5 + <parameter name="PointSpreadFunction" file="$(PACKAGE_SOURCE)/inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h/irf_file.fits"/>
  6 + <parameter name="EnergyDispersion" file="$(PACKAGE_SOURCE)/inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h/irf_file.fits"/>
  7 + <parameter name="Background" file="$(PACKAGE_SOURCE)/inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h/irf_file.fits"/>
7 8 </observation>
8 9 </observation_list>
... ...
inst/cta/test/data/irf_unbinned.xml
  View file @84df48d
1 1 <observation_list title="observation library">
2 2 <observation name="Crab" id="000001" instrument="CTA">
3 3 <parameter name="EventList" file="$(PACKAGE_SOURCE)/inst/cta/test/data/crab_events.fits"/>
4   - <parameter name="EffectiveArea" file="$(PACKAGE_SOURCE)/inst/cta/caldb/dummy/bcf/hard/000001/ea_test.fits"/>
5   - <parameter name="PointSpreadFunction" file="$(PACKAGE_SOURCE)/inst/cta/caldb/dummy/bcf/hard/000001/psf_test.fits"/>
6   - <parameter name="EnergyDispersion" file="$(PACKAGE_SOURCE)/inst/cta/caldb/dummy/bcf/hard/000001/edisp_test.fits"/>
  4 + <parameter name="EffectiveArea" file="$(PACKAGE_SOURCE)/inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h/irf_file.fits"/>
  5 + <parameter name="PointSpreadFunction" file="$(PACKAGE_SOURCE)/inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h/irf_file.fits"/>
  6 + <parameter name="EnergyDispersion" file="$(PACKAGE_SOURCE)/inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h/irf_file.fits"/>
  7 + <parameter name="Background" file="$(PACKAGE_SOURCE)/inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h/irf_file.fits"/>
7 8 </observation>
8 9 </observation_list>
... ...
inst/cta/test/data/obs_binned.xml
  View file @84df48d
... ... @@ -4,5 +4,6 @@
4 4 <parameter name="EffectiveArea" file=""/>
5 5 <parameter name="PointSpreadFunction" file=""/>
6 6 <parameter name="EnergyDispersion" file=""/>
  7 + <parameter name="Background" file=""/>
7 8 </observation>
8 9 </observation_list>
... ...
inst/cta/test/data/obs_unbinned.xml
  View file @84df48d
... ... @@ -4,5 +4,6 @@
4 4 <parameter name="EffectiveArea" file=""/>
5 5 <parameter name="PointSpreadFunction" file=""/>
6 6 <parameter name="EnergyDispersion" file=""/>
  7 + <parameter name="Background" file=""/>
7 8 </observation>
8 9 </observation_list>
... ...