From df53f488232c9d156bd6ba97e3e6e6e999ca6d0c Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Wed, 22 Jun 2016 12:36:43 +0200 Subject: [PATCH] Correct computation of alpha in On/Off observations. --- inst/cta/src/GCTAOnOffObservation.cpp | 62 ++++++++++++++++++++++++++++++-------------------------------- inst/cta/test/test_CTA.py | 18 +++++++++--------- 2 files changed, 39 insertions(+), 41 deletions(-) diff --git a/inst/cta/src/GCTAOnOffObservation.cpp b/inst/cta/src/GCTAOnOffObservation.cpp index 1b38fd8..036bc1e 100644 --- a/inst/cta/src/GCTAOnOffObservation.cpp +++ b/inst/cta/src/GCTAOnOffObservation.cpp @@ -536,7 +536,9 @@ void GCTAOnOffObservation::write(GXmlElement& xml) const * @exception GException::invalid_value * No CTA event list found in CTA observation. * - * @todo Set the PHA spectra to zero before filling. + * Fills events from an observation into the On and Off spectra and cumulates + * their livetime. This method adds to any existing content in the On and + * Off spectra. ***************************************************************************/ void GCTAOnOffObservation::fill(const GCTAObservation& obs) { @@ -1027,14 +1029,20 @@ void GCTAOnOffObservation::free_members(void) /***********************************************************************//** - * @brief Compute alpha parameter (ratio of On to Off effective areas) + * @brief Compute vector of alpha parameters * * @param[in] obs CTA observation. * * @exception GException::invalid_value * No CTA response found in CTA observation. * - * @todo What about the livetime? What about the solid angle of the regions? + * Compute the alpha parameters for all energy bins. The alpha parameter + * gives the ratio between the On and Off region background acceptance + * multiplied by the ratio between On and Off region solid angles. + * + * @todo Implement for general sky regions, not only circles + * @todo Integrate over the sky region instead of taking the values at the + * sky region centre ***************************************************************************/ void GCTAOnOffObservation::compute_alpha(const GCTAObservation& obs) { @@ -1044,7 +1052,7 @@ void GCTAOnOffObservation::compute_alpha(const GCTAObservation& obs) // Continue only if there are spectral bins int nreco = ereco.size(); if (nreco > 0) { - + // Initialise On/Off exposure ratios m_alpha.assign(nreco, 0.0); @@ -1061,8 +1069,6 @@ void GCTAOnOffObservation::compute_alpha(const GCTAObservation& obs) // Get CTA observation pointing direction, zenith, and azimuth GCTAPointing obspnt = obs.pointing(); GSkyDir obsdir = obspnt.dir(); - double zenith = obspnt.zenith(); - double azimuth = obspnt.azimuth(); // Loop over reconstructed energies for (int i = 0; i < nreco; ++i) { @@ -1071,8 +1077,8 @@ void GCTAOnOffObservation::compute_alpha(const GCTAObservation& obs) double logEreco = ereco.elogmean(i).log10TeV(); // Initialise effective area totals - double aon = 0.0; - double aoff = 0.0; + double aon(0.0); + double aoff(0.0); // Loop over ON sky regions for (int m = 0; m < m_on_regions.size(); ++m) { @@ -1085,18 +1091,14 @@ void GCTAOnOffObservation::compute_alpha(const GCTAObservation& obs) // Get centre of circular region GSkyDir centreg = skyreg->centre(); - // Compute position of region centre in instrument - // coordinates (should be replaced by proper method to - // yield theta and phi) - double theta = obsdir.dist(centreg); - double phi = obsdir.posang(centreg); - - // Add up effective area - aon += response->aeff(theta, - phi, - zenith, - azimuth, - logEreco); + // Compute instrument direction from region centre + GCTAInstDir instdir(obspnt.instdir(centreg)); + + // Add up acceptance + aon += (*response->background())(logEreco, + instdir.detx(), + instdir.dety()) * + skyreg->solidangle(); } // endif: sky region is of type GSkyRegionCircle @@ -1113,18 +1115,14 @@ void GCTAOnOffObservation::compute_alpha(const GCTAObservation& obs) // Get centre of circular region GSkyDir centreg = skyreg->centre(); - // Compute position of region centre in instrument - // coordinates (should be replaced by proper method to - // yield theta and phi) - double theta = obsdir.dist(centreg); - double phi = obsdir.posang(centreg); - - // Add up effective area - aoff += response->aeff(theta, - phi, - zenith, - azimuth, - logEreco); + // Compute instrument direction from region centre + GCTAInstDir instdir(obspnt.instdir(centreg)); + + // Add up acceptance + aoff += (*response->background())(logEreco, + instdir.detx(), + instdir.dety()) * + skyreg->solidangle(); } // endif: sky region is of type GSkyRegionCircle diff --git a/inst/cta/test/test_CTA.py b/inst/cta/test/test_CTA.py index 5a3e0c5..233db87 100755 --- a/inst/cta/test/test_CTA.py +++ b/inst/cta/test/test_CTA.py @@ -252,7 +252,7 @@ class Test(gammalib.GPythonTestSuite): crab = gammalib.GSkyDir() crab.radec_deg(83.6331, 22.0145) on = gammalib.GSkyRegions() - on.append(gammalib.GSkyRegionCircle(crab, 0.5)) + on.append(gammalib.GSkyRegionCircle(crab, 0.2)) # Create Off region crab.radec_deg(83.6331, 23.5145) @@ -291,21 +291,21 @@ class Test(gammalib.GPythonTestSuite): # Test On/Off model fitting results sky = outobs.models()['Crab'] bgd = outobs.models()['Background'] - self.test_value(sky['Prefactor'].value(), 6.15281e-16, 1.0e-20, + self.test_value(sky['Prefactor'].value(), 5.95777e-16, 1.0e-20, 'Check sky model prefactor value') - self.test_value(sky['Prefactor'].error(), 2.35727e-17, 1.0e-20, + self.test_value(sky['Prefactor'].error(), 2.02034e-17, 1.0e-20, 'Check sky model prefactor error') - self.test_value(sky['Index'].value(), -2.51599, 1.0e-4, + self.test_value(sky['Index'].value(), -2.50932, 1.0e-4, 'Check sky model index value') - self.test_value(sky['Index'].error(), 0.037628, 1.0e-4, + self.test_value(sky['Index'].error(), 0.0304839, 1.0e-4, 'Check sky model index error') - self.test_value(bgd['Prefactor'].value(), 1.04038, 1.0e-4, + self.test_value(bgd['Prefactor'].value(), 1.09509, 1.0e-4, 'Check background model prefactor value') - self.test_value(bgd['Prefactor'].error(), 0.130852, 1.0e-4, + self.test_value(bgd['Prefactor'].error(), 0.140468, 1.0e-4, 'Check background model prefactor error') - self.test_value(bgd['Index'].value(), 0.539376, 1.0e-4, + self.test_value(bgd['Index'].value(), 0.54909, 1.0e-4, 'Check background model index value') - self.test_value(bgd['Index'].error(), 0.0855112, 1.0e-4, + self.test_value(bgd['Index'].error(), 0.0872644, 1.0e-4, 'Check background model index error') # Save PHA, ARF and RMFs -- libgit2 0.21.4