Use additional ARF columns to store background and alpha vectors

Jürgen Knödlseder
1 parent e8e73b77
Showing 3 changed files with 41 additions and 48 deletions
inst/cta/include/GCTAOnOffObservation.hpp
inst/cta/pyext/GCTAOnOffObservation.i
inst/cta/src/GCTAOnOffObservation.cpp
@@ -88,7 +88,6 @@ public:
     const GPha& on_spec(void) const;
     const GPha& off_spec(void) const;
     const GArf& arf(void) const;
-	const GArf& bgd(void) const;
     const GRmf& rmf(void) const;
     void        fill(const GCTAObservation& obs);
     void        compute_response(const GCTAObservation& obs,
@@ -107,19 +106,17 @@ protected:
 	double model_off(const GModels& models, int ibin, GVector* mod_grad) const;
     // Protected data members
-    std::string         m_instrument;   //!< Instrument name
-    GCTAResponse*       m_response;     //!< Pointer to instrument response functions
-	double              m_ontime;       //!< Ontime (seconds)
-    double              m_livetime;     //!< Livetime (seconds)
-    double              m_deadc;        //!< Deadtime correction (livetime/ontime)
-    GPha 		        m_on_spec;      //!< On counts spectrum
-    GPha 		        m_off_spec;     //!< Off counts spectrum
-    GArf                m_arf;          //!< Effective area vector
-	GArf                m_bgd;          //!< Background rate vector
-    GRmf                m_rmf;          //!< Energy dispersion matrix
-    GSkyRegions         m_on_regions;   //!< Container of On region
-    GSkyRegions         m_off_regions;  //!< Container of Off regions
-    std::vector<double> m_alpha;        //!< Ratios of On/Off exposure
+    std::string         m_instrument;  //!< Instrument name
+    GCTAResponse*       m_response;    //!< Pointer to IRFs
+	double              m_ontime;      //!< Ontime (seconds)
+    double              m_livetime;    //!< Livetime (seconds)
+    double              m_deadc;       //!< Deadtime correction (livetime/ontime)
+    GPha 		        m_on_spec;     //!< On counts spectrum
+    GPha 		        m_off_spec;    //!< Off counts spectrum
+    GArf                m_arf;         //!< Auxiliary Response Function vector
+    GRmf                m_rmf;         //!< Redistribution matrix
+    GSkyRegions         m_on_regions;  //!< Container of On region
+    GSkyRegions         m_off_regions; //!< Container of Off regions
 };
@@ -267,18 +264,6 @@ const GArf&amp; GCTAOnOffObservation::arf(void) const
 /***********************************************************************//**
- * @brief Return Background Rate
- *
- * @return Background Rate.
- ***************************************************************************/
-inline
-const GArf& GCTAOnOffObservation::bgd(void) const
-{
-    return (m_bgd);
-}
-
-
-/***********************************************************************//**
  * @brief Return Redistribution Matrix File
  *
  * @return Redistribution Matrix File.
@@ -291,9 +276,9 @@ const GRmf&amp; GCTAOnOffObservation::rmf(void) const
 /***********************************************************************//**
- * @brief Return Redistribution Matrix File
+ * @brief Return number of observed events
  *
- * @return Redistribution Matrix File.
+ * @return Number of observed events.
  ***************************************************************************/
 inline
 int GCTAOnOffObservation::nobserved(void) const
@@ -95,7 +95,6 @@ public:
     const GPha& on_spec(void) const;
     const GPha& off_spec(void) const;
     const GArf& arf(void) const;
-	const GArf& bgd(void) const;
     const GRmf& rmf(void) const;
     void        fill(const GCTAObservation& obs);
     void        compute_response(const GCTAObservation& obs,
@@ -369,7 +369,7 @@ void GCTAOnOffObservation::read(const GXmlElement&amp; xml)
 			// Get filename
 			std::string filename = par->attribute("file");
-			// load on regions
+			// Load on regions
 			m_on_regions.load(filename);
 			// Increase number of parameters
@@ -395,7 +395,7 @@ void GCTAOnOffObservation::read(const GXmlElement&amp; xml)
 			// Get filename
 			std::string filename = par->attribute("file");
-			// load arf
+			// Load arf
 			m_arf.load(filename);
 			// Increase number of parameters
@@ -408,7 +408,7 @@ void GCTAOnOffObservation::read(const GXmlElement&amp; xml)
 			// Get filename
 			std::string filename = par->attribute("file");
-			// load Rmf
+			// Load Rmf
 			m_rmf.load(filename);
 			// Increase number of parameters
@@ -736,6 +736,9 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
 	// Allocate working array
 	GVector colvar(npars);
+
+    // Get reference to alpha parameters
+    const std::vector<double>& alpha = m_arf["ALPHA"];
 	// Check that there is at least one parameter
 	if (npars > 0) {
@@ -759,7 +762,7 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
             double nbgd = model_off(models, i, &bgd_grad);
             // Skip bin if model is too small (avoids -Inf or NaN gradients)
-            double nonpred = ngam + m_alpha[i] * nbgd;
+            double nonpred = ngam + alpha[i] * nbgd;
             if ((nbgd <= minmod) || (nonpred <= minmod)) {
                 #if defined(G_LIKELIHOOD_DEBUG)
                 n_small_model++;
@@ -783,10 +786,10 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
             // Fill derivatives
             double fa         = non/nonpred;
             double fb         = fa/nonpred;
-            double fc         = m_alpha[i] * fb;
-            double fd         = fc * m_alpha[i] + noff/(nbgd*nbgd);
+            double fc         = alpha[i] * fb;
+            double fd         = fc * alpha[i] + noff/(nbgd*nbgd);
             double sky_factor = 1.0 - fa;
-            double bgd_factor = 1.0 + m_alpha[i] - m_alpha[i] * fa - noff/nbgd;
+            double bgd_factor = 1.0 + alpha[i] - alpha[i] * fa - noff/nbgd;
             // Loop over all parameters
             for (int j = 0; j < npars; ++j) {
@@ -977,11 +980,9 @@ void GCTAOnOffObservation::init_members(void)
     m_on_spec.clear();
     m_off_spec.clear();
     m_arf.clear();
-	m_bgd.clear();
     m_rmf.clear();
     m_on_regions.clear();
     m_off_regions.clear();
-    m_alpha.clear();
     // Return
     return;
@@ -1003,11 +1004,9 @@ void GCTAOnOffObservation::copy_members(const GCTAOnOffObservation&amp; obs)
     m_on_spec     = obs.m_on_spec;
     m_off_spec    = obs.m_off_spec;
     m_arf         = obs.m_arf;
-	m_bgd         = obs.m_bgd;
     m_rmf         = obs.m_rmf;
     m_on_regions  = obs.m_on_regions;
     m_off_regions = obs.m_off_regions;
-    m_alpha       = obs.m_alpha;
     // Clone members
     m_response = (obs.m_response != NULL) ? obs.m_response->clone() : NULL;
@@ -1054,7 +1053,7 @@ void GCTAOnOffObservation::compute_alpha(const GCTAObservation&amp; obs)
     if (nreco > 0) {
 		// Initialise On/Off exposure ratios
-        m_alpha.assign(nreco, 0.0);
+        std::vector<double> alpha(nreco, 0.0);
         // Get CTA response pointer. Throw an exception if no response is found
         const GCTAResponseIrf* response =
@@ -1130,10 +1129,13 @@ void GCTAOnOffObservation::compute_alpha(const GCTAObservation&amp; obs)
 			// Compute alpha for this energy bin
 			if (aoff > 0.0) {
-				m_alpha[i] = aon/aoff;
+				alpha[i] = aon/aoff;
 			}
         } // endfor: looped over reconstructed energies
+
+        // Append alpha vector to ARF
+        m_arf.append("ALPHA", alpha);
     } // endif: there were energy bins
@@ -1246,15 +1248,15 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs)
     // Continue only if there are spectral bins
     int nreco = ereco.size();
     if (nreco > 0) {
+
+		// Initialise background rayes
+        std::vector<double> background(nreco, 0.0);
 		// Get CTA observation pointing direction, zenith, and azimuth
 		GCTAPointing obspnt(obs.pointing());
 		double       zenith  = obspnt.zenith();
 		double       azimuth = obspnt.azimuth();
-		// Initialize background rate array
-    	m_bgd = GArf(ereco);
-	
         // Get pointer on CTA IRF response
         const GCTAResponseIrf* rsp =
               dynamic_cast<const GCTAResponseIrf*>(obs.response());
@@ -1296,12 +1298,16 @@ void GCTAOnOffObservation::compute_bgd(const GCTAObservation&amp; obs)
                 double logEreco = m_on_spec.ebounds().elogmean(i).log10TeV();
                 // Get rate in events/s/MeV
-                m_bgd[i] += (*bgd)(logEreco, offdir.detx(), offdir.dety()) *
-                            skyreg->solidangle();
+                background[i] += (*bgd)(logEreco,
+                                        offdir.detx(),
+                                        offdir.dety()) * skyreg->solidangle();
             } // endfor: looped over energy bins
         } // endfor: looped over sky regions
+
+        // Append background vector to ARF
+        m_arf.append("BACKGROUND", background);
 	} // endif: there were spectral bins
@@ -1557,6 +1563,9 @@ double GCTAOnOffObservation::model_off(const GModels&amp; models,
 	// Continue only if bin number is in range and if there are parameters
 	if ((ibin < m_off_spec.size()) && (npars > 0))  {
+
+        // Get reference to background rates
+        const std::vector<double>& background = m_arf["BACKGROUND"];
         // Get energy bin mean and width
         const GEnergy emean  = m_off_spec.ebounds().elogmean(ibin);
@@ -1564,7 +1573,7 @@ double GCTAOnOffObservation::model_off(const GModels&amp; models,
         // Compute normalisation factor
         double exposure = m_off_spec.exposure();
-        double norm     = m_bgd[ibin] * exposure * ewidth;
+        double norm     = background[ibin] * exposure * ewidth;
         // Loop over models
         for (int j = 0; j < models.size(); ++j) {