Merge branch '1371-adapt-ctools-for-new-binned-interface' into integration

Conflicts: ChangeLog NEWS src/ctbkgcube/ctbkgcube.cpp test/test_ctools.sh

Merge branch '1371-adapt-ctools-for-new-binned-interface' into integration
Conflicts: ChangeLog NEWS src/ctbkgcube/ctbkgcube.cpp test/test_ctools.sh
Jürgen Knödlseder
2 parents 0cb347f7 42aae297
Showing 23 changed files with 348 additions and 529 deletions
ChangeLog
NEWS
doc/source/reference_manual/cspull.rst
doc/source/reference_manual/ctbutterfly.rst
doc/source/reference_manual/ctlike.rst
doc/source/reference_manual/ctmodel.rst
doc/source/reference_manual/cttsmap.rst
doc/source/reference_manual/ctulimit.rst
pyext/ctbkgcube.i
scripts/cspull.py
src/ctbkgcube/ctbkgcube.cpp
src/ctbkgcube/ctbkgcube.hpp
src/ctbkgcube/ctbkgcube.par
src/ctbutterfly/ctbutterfly.par
src/ctlike/ctlike.par
src/ctmodel/ctmodel.par
src/cttsmap/cttsmap.par
src/ctulimit/ctulimit.par
src/support/ctool.cpp
test/data/model_background_cube_const.xml
@@ -11,6 +11,7 @@
         Add ctulimit tool (#1363)
         Allow ctmodel usage for stacked analysis from command line
         Correct rate check in ctobssim for map cubes (#1452)
+        Handle background cubes for stacked analysis (#1371)
 2015-01-08  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
@@ -2,7 +2,7 @@ New Features and Important Changes in ctools 1.0.0
 Author(s) : Juergen Knoedlseder
-23 March 2015
+28 March 2015
 1. Introduction
@@ -50,12 +50,14 @@ None
 8. ctlike - Maximum likelihood fitting
 --------------------------------------
-None
+The bkgcube parameter has been added to allow specification of a background
+cube for stacked analysis.
 9. cttsmap - Generation of Test Statistics map
 ----------------------------------------------
-None
+The bkgcube parameter has been added to allow specification of a background
+cube for stacked analysis.
 10. ctskymap - CTA sky mapping
@@ -71,6 +73,9 @@ the original IRF information. This makes generation of model cubes substantially
 faster. The old method of accessing the IRF response for model computation is
 still possible.
+The bkgcube parameter has been added to allow specification of a background
+cube for stacked analysis.
+
 12. ctexpcube - CTA exposure cube generation
 --------------------------------------------
@@ -94,13 +99,17 @@ None
 16. ctbutterfly - Generate butterfly diagrams
 ---------------------------------------------
-None
+The bkgcube parameter has been added to allow specification of a background
+cube for stacked analysis.
 17. ctulimit - Compute upper limit
 ----------------------------------
 Added an initial version of the upper limit computation tool.
+The bkgcube parameter has been added to allow specification of a background
+cube for stacked analysis.
+
 18. Examples
 ------------
@@ -26,9 +26,15 @@ General parameters
 ``outfile [file]``
     ASCII file containing the individual pull values.
-``ntrials [integer]``
-    Number of samples for generating the pull distribution.
- 	 	 
+``expcube [file]``
+    Exposure cube file (only needed for stacked analysis).
+
+``psfcube [file]``
+    PSF cube file (only needed for stacked analysis).
+
+``bkgcube [file]``
+    Background cube file (only needed for stacked analysis).
+
 ``caldb [string]``
     Calibration database.
@@ -38,6 +44,9 @@ General parameters
 ``(edisp = no) [boolean]``
     Apply energy dispersion to response computation.
+``ntrials [integer]``
+    Number of samples for generating the pull distribution.
+ 	 	 
 ``ra [real]``
     Right Ascension of CTA pointing (J2000, in degrees).
@@ -25,12 +25,15 @@ General parameters
 ``srcname [string]``
     Name of model component for which butterfly should be computed.
-``expcube = NONE [file]``
+``expcube [file]``
     Exposure cube file (only needed for stacked analysis).
-``psfcube = NONE [file]``
+``psfcube [file]``
     PSF cube file (only needed for stacked analysis).
+``bkgcube [file]``
+    Background cube file (only needed for stacked analysis).
+
 ``caldb [string]``
     Calibration database.
@@ -34,6 +34,9 @@ General parameters
 ``psfcube [file]``
     PSF cube file (only needed for stacked analysis).
+``bkgcube [file]``
+    Background cube file (only needed for stacked analysis).
+
 ``caldb [string]``
     Calibration database.
@@ -43,6 +43,9 @@ General parameters
 ``psfcube [file]``
     PSF cube file (only needed for stacked analysis).
+``bkgcube [file]``
+    Background cube file (only needed for stacked analysis).
+
 ``caldb [string]``
     Calibration database.
@@ -31,6 +31,9 @@ General parameters
 ``psfcube [file]``
     PSF cube file (only needed for stacked analysis).
+``bkgcube [file]``
+    Background cube file (only needed for stacked analysis).
+
 ``caldb [string]``
     Calibration database.
@@ -24,12 +24,15 @@ General parameters
 ``srcname [string]``
     Name of model component for which upper limit should be computed.
-``expcube = NONE [file]``
+``expcube [file]``
     Exposure cube file (only needed for stacked analysis).
-``psfcube = NONE [file]``
+``psfcube [file]``
     PSF cube file (only needed for stacked analysis).
+``bkgcube [file]``
+    Background cube file (only needed for stacked analysis).
+
 ``caldb [string]``
     Calibration database.
 /***************************************************************************
  *               ctbkgcube - Background cube generation tool               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2014 by Chia-Chun Lu                                     *
+ *  copyright (C) 2014-2015 by Chia-Chun Lu                                *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -44,11 +44,11 @@ public:
     virtual ~ctbkgcube(void);
     // Methods
-    void                 clear(void);
-    void                 run(void);
-    void                 save(void);
-    const GCTAEventCube& cube(void) const;
-    const GModels&       models(void) const;
+    void                      clear(void);
+    void                      run(void);
+    void                      save(void);
+    const GCTACubeBackground& cube(void) const;
+    const GModels&            models(void) const;
 };
 /***********************************************************************//**
@@ -98,8 +98,11 @@ class cspull(ctools.cscript):
             pars.append(gammalib.GApplicationPar("outfile","f","a","pull.dat","","","Output file name"))
             pars.append(gammalib.GApplicationPar("expcube","s","a","NONE","","","Exposure cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("psfcube","s","a","NONE","","","PSF cube file (only needed for stacked analysis)"))
+            pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Background cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
             pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
+            pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
+            pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials"))
             pars.append(gammalib.GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)"))
             pars.append(gammalib.GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)"))
             pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit (TeV)"))
@@ -107,14 +110,12 @@ class cspull(ctools.cscript):
             pars.append(gammalib.GApplicationPar("tmin","r","h","0.0","","","Start time (MET in s)"))
             pars.append(gammalib.GApplicationPar("tmax","r","a","1800.0","","","Duration (in s)"))
             pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
-            pars.append(gammalib.GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
-            pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
+            pars.append(gammalib.GApplicationPar("npix","i","a","200","","","Number of pixels for binned"))
+            pars.append(gammalib.GApplicationPar("binsz","r","a","0.05","","","Pixel size for binned (deg/pixel)"))
             pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","0","1","Deadtime correction factor"))
             pars.append(gammalib.GApplicationPar("rad","r","h","5.0","0","180","Radius of ROI (deg)"))         
             pars.append(gammalib.GApplicationPar("pattern","s","h","single","","","Observation pattern (single/four)"))
             pars.append(gammalib.GApplicationPar("offset","r","h","1.5","","","Observation pattern offset (deg)"))
-            pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials"))
-            pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
             pars.append_standard()
             pars.append(gammalib.GApplicationPar("logfile","f","h","cspull.log","","","Log filename"))
             pars.save(parfile)
@@ -140,7 +141,7 @@ class cspull(ctools.cscript):
         # Get number of energy bins
         self.m_enumbins = self["enumbins"].integer()
-        
+   
         # Read parameters for binned if requested
         if not self.m_enumbins == 0:
             self.m_npix  = self["npix"].integer()
@@ -34,7 +34,6 @@
 /* __ Method name definitions ____________________________________________ */
 #define G_GET_PARAMETERS                        "ctbkgcube::get_parameters()"
-#define G_FILL_CUBE                  "ctbkgcube::fill_cube(GCTAObservation*)"
 /* __ Debug definitions __________________________________________________ */
@@ -209,9 +208,6 @@ void ctbkgcube::run(void)
         log.cout(true);
     }
-    // Initialise run parameters
-    m_livetime = 0.0;
-
     // Get task parameters
     get_parameters();
@@ -248,15 +244,14 @@ void ctbkgcube::run(void)
     // Copy models from observation container and reset output model
     // container
-    m_bkgmdl = m_obs.models();
+    GModels models_orig = m_obs.models();
+    m_bkgmdl            = m_obs.models();
     m_outmdl.clear();
-    m_cube_model = -1;
     // Remove all models that are not CTA background models from the
     // container and put all removed components in the output
     // container
     int num = m_bkgmdl.size();
-
     for (int i = num-1; i >= 0; --i) {
         // Flag removal
@@ -309,75 +304,28 @@ void ctbkgcube::run(void)
         log.header1("Generate background cube");
     }
-    // Initialise observation counter
-    int n_observations = 0;
-
-    // Loop over all observations in the container
-    for (int i = 0; i < m_obs.size(); ++i) {
-      
-        // Get CTA observation
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-
-        // Continue only if observation is a CTA observation
-        if (obs != NULL) {
-
-            // Write header for observation
-            if (logTerse()) {
-                if (obs->name().length() > 1) {
-                    log.header3("Observation "+obs->name());
-                }
-                else {
-                    log.header3("Observation "+gammalib::str(i));
-                }
-            }
-
-            // Increment number of observations
-            n_observations++;
-
-            // Fill the cube
-            fill_cube(obs);
-
-        } // endif: CTA observation found
+    // Assign background models to container
+    m_obs.models(m_bkgmdl);
-    } // endfor: looped over observations
+    // Fill background cube from observations
+    m_background.fill(m_obs);
-    // Log results
-    if (logTerse()) {
-        log << gammalib::parformat("Accumulated livetime");
-        log << m_livetime << " sec" << std::endl;
-    }
+    // Create a background model for the output background cube and append
+    // that model to the input model in place of the original
+    // background models
+    // TODO: We might think of creating the spectral model via user parameter
+    GModelSpectralPlaw      spectral(1.0, 0.0, GEnergy(1.0, "TeV"));
+    GCTAModelCubeBackground model(spectral);
-    // Re-normalize cube to get units of counts/MeV/s/sr
-    if (m_livetime > 0.0) {
-    
-        // Loop over all bins in background cube and divide the content
-        // by the total livetime.
-        for (int i = 0; i < m_bkgcube.size(); ++i) {
-            GCTAEventBin* bin  = m_bkgcube[i];
-            double        rate = bin->counts() / m_livetime;
-            bin->counts(rate);
-        }
+    // Set model name
+    model.name("BackgroundModel");
-    } // endif: livetime was positive
+    // Set model instrument
+    // TODO: Account for possibility to have observations from different
+    // IACTs in the same container
+    model.instruments("CTA,HESS,MAGIC,VERITAS");
-    // Log results
-    if (logTerse()) {
-        log << gammalib::parformat("Number of observations");
-        log << n_observations << std::endl;
-    }
-
-    // Create a background model for output background cube and append
-    // that model to the input model in place of the original
-    // background model
-    GEnergies energies;
-    for (int i = 0; i < m_bkgcube.ebins(); ++i) {
-        energies.append(m_bkgcube.energy(i));
-    }
-    GModelSpatialDiffuseCube spatial(m_bkgcube.map(), energies, 1.0);
-    GModelSpectralPlaw       spectral(1.0, 0.0, GEnergy(1.0, "TeV"));
-    GCTAModelCubeBackground  model(spatial, spectral);
-    model.name("ctbkgcube default background model");
-    m_cube_model = m_outmdl.size(); // Store the slot number for save()
+    // Append model to output container
     m_outmdl.append(model);
     // Log output model
@@ -387,6 +335,9 @@ void ctbkgcube::run(void)
         log << m_outmdl << std::endl;
     }
+    // Recover original models
+    m_obs.models(models_orig);
+
     // Return
     return;
 }
@@ -405,47 +356,19 @@ void ctbkgcube::save(void)
         log.header1("Save background cube");
     }
-    // Create energies container from energy boundaries
-    GEnergies energies;
-    for (int i = 0; i < m_ebounds.size(); ++i) {
-        energies.append(m_ebounds.elogmean(i));
-    }
-
     // Get output filenames
     std::string outfile  = (*this)["outcube"].filename();
     std::string outmodel = (*this)["outmodel"].filename();
-    // Create empty FITS file
-    GFits fits;
-
-    // Write background cube
-    m_bkgcube.map().write(fits);
-
-    // Write energies
-    energies.write(fits);
-    
-    // Save FITS file
-    fits.saveto(outfile, clobber());
+    // Save background cube
+    m_background.save(outfile, clobber());
     // Write output models if filename is valid
     if ((outmodel.length() > 0) && (gammalib::tolower(outmodel) != "none")) {
-        // Set filename of map cube.
-        if (m_cube_model >= 0) {
-            GCTAModelCubeBackground* model =
-                dynamic_cast<GCTAModelCubeBackground*>(m_outmdl[m_cube_model]);
-            if (model != NULL) {
-                GModelSpatialDiffuseCube* spatial =
-                    dynamic_cast<GModelSpatialDiffuseCube*>(model->spatial());
-                if (spatial != NULL) {
-                    spatial->filename(outfile);
-                }
-            }
-        }
-
-        // Save output model container
+        // Save output model for binned analyses
         m_outmdl.save(outmodel);
-        
+
     }
     // Return
@@ -468,12 +391,9 @@ void ctbkgcube::init_members(void)
     m_outcube.clear();
     m_outmodel.clear();
     m_obs.clear();
-    m_bkgcube.clear();
+    m_background.clear();
     m_bkgmdl.clear();
     m_outmdl.clear();
-    m_ebounds.clear();
-    m_livetime   = 0.0;
-    m_cube_model = -1;
     // Return
     return;
@@ -491,12 +411,9 @@ void ctbkgcube::copy_members(const ctbkgcube&amp; app)
     m_outmodel   = app.m_outmodel;
     m_outcube    = app.m_outcube;
     m_obs        = app.m_obs;
-    m_bkgcube    = app.m_bkgcube;
+    m_background = app.m_background;
     m_bkgmdl     = app.m_bkgmdl;
     m_outmdl     = app.m_outmdl;
-    m_ebounds    = app.m_ebounds;
-    m_livetime   = app.m_livetime;
-    m_cube_model = app.m_cube_model;
     // Return
     return;
@@ -536,22 +453,6 @@ void ctbkgcube::get_parameters(void)
     } // endif: there was no observation in the container
-
-    // If there are no models associated with the observations then load now
-    // the model definition from the XML file
-    if (m_obs.models().size() == 0) {
-        if ((*this)["inmodel"].is_undefined()) {
-            std::string msg = "No model definition XML file specified. "
-                              "Please set the \"inmodel\" parameter to the "
-                              "XML file that contains the background model "
-                              "definition.";
-            throw GException::invalid_value(G_GET_PARAMETERS, msg);
-        }
-        std::string inmodel = (*this)["inmodel"].filename();
-        GModels     models(inmodel);
-        m_obs.models(models);
-    }
-
     // Get the incube filename
     std::string incube = (*this)["incube"].filename();
@@ -560,28 +461,35 @@ void ctbkgcube::get_parameters(void)
         (gammalib::strip_whitespace(incube) == "")) {
         // Create an event cube based on task parameters
-        m_bkgcube = create_cube(m_obs);
+        GCTAEventCube cube = create_cube(m_obs);
+
+        // Define background cube
+        m_background = GCTACubeBackground(cube);
     } // endif: filename was not valid
     // ... otherwise setup the background cube from the counts map
     else {
-    
-        // Load event cube from filename
-        m_bkgcube.load(incube);
-        // Loop over all bins in background cube
-        for (int i = 0; i < m_bkgcube.size(); ++i) {
-            // Get event bin
-            GCTAEventBin* bin = m_bkgcube[i];
-            bin->counts(0.0);
-
-        }
+        // Define background cube from file
+        m_background = GCTACubeBackground(incube);
     } // endelse: cube was loaded from file
-    // Get energy definition
-    m_ebounds = m_bkgcube.ebounds();
+    // If there are no models associated with the observations then load now
+    // the model definition from the XML file
+    if (m_obs.models().size() == 0) {
+        if ((*this)["inmodel"].is_undefined()) {
+            std::string msg = "No model definition XML file specified. "
+                              "Please set the \"inmodel\" parameter to the "
+                              "XML file that contains the background model "
+                              "definition.";
+            throw GException::invalid_value(G_GET_PARAMETERS, msg);
+        }
+        std::string inmodel = (*this)["inmodel"].filename();
+        GModels     models(inmodel);
+        m_obs.models(models);
+    }
     // Read output filenames (if needed)
     if (read_ahead()) {
@@ -592,92 +500,3 @@ void ctbkgcube::get_parameters(void)
     // Return
     return;
 }
-
-/***********************************************************************//**
- * @brief Generate background cube
- *
- * @param[in] obs Pointer to CTA observation.
- *
- * Fills the background cube with the model value for a CTA observation.
- * This method will not need to load the event file as it gets the event
- * ROI and GTI directly from the observation.
- ***************************************************************************/
-void ctbkgcube::fill_cube(GCTAObservation* obs)
-{
-
-    // Continue only if observation pointer is valid
-    if (obs != NULL) {
-
-        // Initialise statistics
-        double sum            = 0.0;
-        int    n_bins_inside  = 0;
-        int    n_bins_outside = 0;
-
-        // Extract region of interest from CTA observation
-        GCTARoi roi = obs->roi();
-
-        // Set GTI of actual observations as the GTI of the event cube
-        m_bkgcube.gti(obs->gti());
-
-        // Get observation livetime
-        double livetime = obs->livetime();
-
-        // Loop over all bins in background cube
-        for (int i = 0; i < m_bkgcube.size(); ++i) {
-
-            // Get event bin
-            GCTAEventBin* bin = m_bkgcube[i];
-
-            // Continue only if binned in contained in ROI
-            if (roi.contains(*bin)) {
-            
-                // Compute model value for event bin. The model value is
-                // given in counts/MeV/s/sr.
-                double model = m_bkgmdl.eval(*bin, *obs);
-
-                // Compute number events
-                sum += model * bin->size();
-
-                // Multiply by livetime to get the correct weighting for
-                // each observation. We divide by the total livetime later
-                // to get the background model in units of counts/MeV/s/sr.
-                model *= livetime;
-
-                // Add existing number of counts
-                model += bin->counts();
-
-                // Store cumulated value (units: counts/MeV/sr)
-                bin->counts(model);
-
-                // Increment contained bin counter
-                n_bins_inside++;
-                
-            } // endif: bin was contained in RoI
-            
-            // ... otherwise count bin as an outsider
-            else {
-                n_bins_outside++;
-            }
-
-        } // endfor: looped over all bins
-
-        // Accumulate livetime
-        m_livetime += livetime;
-
-        // Log results
-        if (logTerse()) {
-            log << gammalib::parformat("Bins within RoI");
-            log << n_bins_inside << std::endl;
-            log << gammalib::parformat("Bins outside RoI");
-            log << n_bins_outside << std::endl;
-            log << gammalib::parformat("Background events in cube");
-            log << sum << std::endl;
-            log << gammalib::parformat("Cube livetime");
-            log << livetime << " sec" << std::endl;
-        }
-
-    } // endif: observation pointer was not valid
-
-    // Return
-    return;
-}
@@ -56,11 +56,11 @@ public:
     ctbkgcube& operator=(const ctbkgcube& app);
     // Methods
-    void                 clear(void);
-    void                 run(void);
-    void                 save(void);
-    const GCTAEventCube& cube(void) const;
-    const GModels&       models(void) const;
+    void                      clear(void);
+    void                      run(void);
+    void                      save(void);
+    const GCTACubeBackground& cube(void) const;
+    const GModels&            models(void) const;
 protected:
     // Protected methods
@@ -71,31 +71,29 @@ protected:
     void fill_cube(GCTAObservation* obs);
     // Protected members
-    std::string   m_outcube;     //!< Filename of output cube
-    std::string   m_outmodel;    //!< Filename of output XML model
-    GObservations m_obs;         //!< Observation container
-    GCTAEventCube m_bkgcube;     //!< Background cube
-    GModels       m_bkgmdl;      //!< CTA background models
-    GModels       m_outmdl;      //!< Output models
-    GEbounds      m_ebounds;     //!< Energy boundaries
-    double        m_livetime;    //!< Total livetime for normalization
-    int           m_cube_model;  //!< Cube model index
+    std::string        m_outcube;     //!< Filename of output cube
+    std::string        m_outmodel;    //!< Filename of output XML model
+    GObservations      m_obs;         //!< Observation container
+    GCTACubeBackground m_background;  //!< Background cube response
+    GModels            m_bkgmdl;      //!< CTA background models
+    GModels            m_outmdl;      //!< Output models
 };
 /***********************************************************************//**
- * @brief Return event cube containing background rate
+ * @brief Return background response cube containing background rate
  *
- * @return Background event cube.
+ * @return Background response cube.
  ***************************************************************************/
 inline
-const GCTAEventCube& ctbkgcube::cube(void) const
+const GCTACubeBackground& ctbkgcube::cube(void) const
 {
     // Return background cube
-    return (m_bkgcube);
+    return (m_background);
 }
+
 /***********************************************************************//**
  * @brief Return background model container
  *
@@ -18,10 +18,10 @@
 # Observational parameters
 #===================
 inobs,    f, a, "NONE",,, "Input event list or observation definition file"
-inmodel,  f, a, "NONE",,, "Input (background) model XML file"
+inmodel,  f, a, "NONE",,, "Input model XML file"
 incube,   f, a, "NONE",,, "Counts cube for background cube definition"
 outcube,  f, a, "bkgcube.fits",,, "Output background cube file"
-outmodel, f, a, "NONE",,, "Output (background) model XML file"
+outmodel, f, a, "NONE",,, "Output model XML file"
 caldb,    s, a, "dummy",,, "Calibration database"
 irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
 #
@@ -22,6 +22,7 @@ inmodel,  f, a, &quot;$CTOOLS/share/models/crab.xml&quot;,,, &quot;Source model&quot;
 srcname,  s, a, "Crab",,,"Source of interest"
 expcube,  f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
 psfcube,  f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
+bkgcube,  f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
 caldb,    s, a, "dummy",,, "Calibration database"
 irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
 outfile,  f, a, "butterfly.txt",,, "Output ascii file"
@@ -21,6 +21,7 @@ inobs,    f, a, &quot;events.fits&quot;,,, &quot;Event list, counts cube or observation definit
 inmodel,  f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
 expcube,  f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
 psfcube,  f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
+bkgcube,  f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
 caldb,    s, a, "dummy",,, "Calibration database"
 irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
 outmodel, f, a, "crab_results.xml",,, "Source model output file"
@@ -27,6 +27,7 @@ edisp,    b, h, no,,, &quot;Apply energy dispersion?&quot;
 inobs,    f, a, "NONE",,, "Event list, counts cube or observation definition file"
 expcube,  f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
 psfcube,  f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
+bkgcube,  f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
 caldb,    s, a, "dummy",,, "Calibration database"
 irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
 ra,       r, a, 83.63,0,360, "RA of pointing (degrees)"
@@ -22,6 +22,7 @@ inmodel, f, a, &quot;$CTOOLS/share/models/crab.xml&quot;,,, &quot;Source model&quot;
 srcname, s, a, "Crab",,, "Test source"
 expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
 psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
+bkgcube, f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
 caldb,   s, a, "dummy",,, "Calibration database"
 irf,     s, a, "cta_dummy_irf",,, "Instrument response function"
 outmap,  f, a, "tsmap.fits",,, "Output Test Statistic map"
@@ -22,6 +22,7 @@ inmodel,  f, a, &quot;$CTOOLS/share/models/crab.xml&quot;,,, &quot;Source model&quot;
 srcname,  s, a, "Crab",,, "Source of interest"
 expcube,  f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
 psfcube,  f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
+bkgcube,  f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
 caldb,    s, a, "dummy",,, "Calibration database"
 irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
@@ -732,23 +732,26 @@ void ctool::set_obs_response(GCTAObservation* obs)
     // If the observation contains a counts cube, then first check whether
     // the expcube and psfcube parameters exist and are not NONE
     if (dynamic_cast<const GCTAEventCube*>(obs->events()) != NULL) {
-        if (has_par("expcube") && has_par("psfcube")) {
+        if (has_par("expcube") && has_par("psfcube") && has_par("bkgcube")) {
             // Get filenames
             std::string expcube = (*this)["expcube"].filename();
             std::string psfcube = (*this)["psfcube"].filename();
+            std::string bkgcube = (*this)["bkgcube"].filename();
             // Extract response information if available
-            if ((expcube != "NONE") && (psfcube != "NONE") &&
+            if ((expcube != "NONE") && (psfcube != "NONE") && (bkgcube != "NONE") &&
                 (gammalib::strip_whitespace(expcube) != "") &&
-                (gammalib::strip_whitespace(psfcube) != "")) {
+                (gammalib::strip_whitespace(psfcube) != "") &&
+                (gammalib::strip_whitespace(bkgcube) != "")) {
                 // Get exposure and PSF cubes
-                GCTACubeExposure exposure(expcube);
-                GCTACubePsf      psf(psfcube);
+                GCTACubeExposure   exposure(expcube);
+                GCTACubePsf        psf(psfcube);
+                GCTACubeBackground background(bkgcube);
                 // Set reponse
-                obs->response(exposure, psf);
+                obs->response(exposure, psf, background);
                 // Signal response availability
                 has_response = true;
-<?xml version="1.0" standalone="no"?>	
-<source_library title="source library">
-  <source name="CTABackgroundModel" type="CTACubeBackground" instrument="CTA">
-    <spectrum type="ConstantValue">
-      <parameter name="Value" scale="1.0" min="0.01" max="1000.0" value="1.0" free="1"/>	
-    </spectrum>
-    <spatialModel type="MapCubeFunction" file="data/test_cube.fits">
-      <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
-    </spatialModel>
-  </source>	
-</source_library>
@@ -6,8 +6,5 @@
       <parameter name="Index"     scale="1.0"  value="0.0"  min="-5.0" max="+5.0"   free="1"/>	
       <parameter name="Scale"     scale="1e6"  value="1.0"  min="0.01" max="1000.0" free="0"/>	
     </spectrum>
-    <spatialModel type="MapCubeFunction" file="data/test_cube.fits">
-      <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
-    </spatialModel>
   </source>	
 </source_library>
@@ -91,27 +91,6 @@ else
   echo " events.fits file is not found"
   exit 1
 fi
-#
-# Run 2
-$ctobssim inmodel="data/model_background_cube.xml" \
-          outevents="events2.fits" \
-          caldb="irf" \
-          irf="cta_dummy_irf" \
-          ra=84.17263 \
-          dec=22.01444 \
-          rad=1.0 \
-          tmin=0.0 \
-          tmax=1.0 \
-          emin=0.3 \
-          emax=8.0 chatter=4
-$ECHO -n "."
-if [ -s "events2.fits" ]
-then
-  $ECHO -n "."
-else
-  echo " events2.fits file is not found"
-  exit 1
-fi
 $ECHO " ok"
@@ -148,7 +127,7 @@ $ECHO -n &quot;Test ctbin: &quot;
 #
 # Run 1
 $ctbin inobs="events.fits" \
-       outcube="cntmap.fits" \
+       outcube="cntmap1.fits" \
        emin=0.1 \
        emax=100.0 \
        enumbins=20 \
@@ -161,27 +140,27 @@ $ctbin inobs=&quot;events.fits&quot; \
        yref=22.01 \
        proj="CAR"
 $ECHO -n "."
-if [ -s "cntmap.fits" ]
+if [ -s "cntmap1.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " cntmap.fits file is not found"
+  $ECHO " cntmap1.fits file is not found"
   exit 1
 fi
 #
 # Run 2
-$ctbin inobs="events2.fits" \
+$ctbin inobs="events.fits" \
        outcube="cntmap2.fits" \
        emin=0.3 \
        emax=8.0 \
-       enumbins=10 \
+       enumbins=7 \
        ebinalg="LOG" \
-       nxpix=50 \
-       nypix=50 \
-       binsz=0.02 \
+       nxpix=30 \
+       nypix=30 \
+       binsz=0.05 \
        coordsys="CEL" \
-       xref=84.17263 \
-       yref=22.01444 \
+       xref=86.63 \
+       yref=22.01 \
        proj="CAR"
 $ECHO -n "."
 if [ -s "cntmap2.fits" ]
@@ -193,6 +172,168 @@ else
 fi
 $ECHO " ok"
+#
+# Test ctexpcube
+# ==============
+$ECHO -n "Test ctexpcube: "
+#
+# Run 1
+$ctexpcube inobs="data/crab_events.fits" \
+           incube="NONE" \
+           outcube="expcube1.fits" \
+           caldb="irf" \
+           irf="cta_dummy_irf" \
+           emin=0.1 \
+           emax=100.0 \
+           enumbins=20 \
+           nxpix=200 \
+           nypix=200 \
+           binsz=0.02 \
+           coordsys="CEL" \
+           xref=83.63 \
+           yref=22.01 \
+           proj="CAR"
+$ECHO -n "."
+if [ -s "expcube1.fits" ]
+then
+$ECHO -n "."
+else
+$ECHO " expcube1.fits file is not found"
+exit 1
+fi
+#
+# Run 2
+$ctexpcube inobs="events.fits" \
+           incube="cntmap2.fits" \
+           outcube="expcube2.fits" \
+           caldb="irf" \
+           irf="cta_dummy_irf"
+$ECHO -n "."
+if [ -s "expcube2.fits" ]
+then
+$ECHO -n "."
+else
+$ECHO " expcube2.fits file is not found"
+exit 1
+fi
+$ECHO " ok"
+
+
+#
+# Test ctpsfcube
+# ==============
+$ECHO -n "Test ctpsfcube: "
+#
+# Run 1
+$ctpsfcube inobs="data/crab_events.fits" \
+           incube="NONE" \
+           outcube="psfcube1.fits" \
+           caldb="irf" \
+           irf="cta_dummy_irf" \
+           emin=0.1 \
+           emax=100.0 \
+           enumbins=20 \
+           nxpix=10 \
+           nypix=10 \
+           binsz=0.4 \
+           coordsys="CEL" \
+           xref=83.63 \
+           yref=22.01 \
+           proj="CAR" \
+           amax=0.3 \
+           anumbins=10
+$ECHO -n "."
+if [ -s "psfcube1.fits" ]
+then
+$ECHO -n "."
+else
+$ECHO " psfcube1.fits file is not found"
+exit 1
+fi
+#
+# Run 2
+$ctpsfcube inobs="events.fits" \
+           incube="cntmap2.fits" \
+           outcube="psfcube2.fits" \
+           caldb="irf" \
+           irf="cta_dummy_irf" \
+           amax=0.3 \
+           anumbins=2
+$ECHO -n "."
+if [ -s "psfcube2.fits" ]
+then
+$ECHO -n "."
+else
+$ECHO " psfcube2.fits file is not found"
+exit 1
+fi
+$ECHO " ok"
+
+
+#
+# Test ctbkgcube
+# ==============
+$ECHO -n "Test ctbkgcube: "
+#
+# Run 1
+$ctbkgcube inobs="data/crab_events.fits" \
+           inmodel="data/crab.xml" \
+           incube="NONE" \
+           outcube="bkgcube1.fits" \
+           outmodel="bkgcube1.xml" \
+           caldb="irf" \
+           irf="cta_dummy_irf" \
+           emin=0.1 \
+           emax=100.0 \
+           enumbins=20 \
+           nxpix=10 \
+           nypix=10 \
+           binsz=0.4 \
+           coordsys="CEL" \
+           xref=83.63 \
+           yref=22.01 \
+           proj="CAR"
+$ECHO -n "."
+if [ -s "bkgcube1.fits" ]
+then
+$ECHO -n "."
+else
+$ECHO " bkgcube1.fits file is not found"
+exit 1
+fi
+if [ -s "bkgcube1.xml" ]
+then
+$ECHO -n "."
+else
+$ECHO " bkgcube1.xml file is not found"
+exit 1
+fi
+#
+# Run 2
+$ctbkgcube inobs="events.fits" \
+           inmodel="data/crab.xml" \
+           incube="cntmap2.fits" \
+           caldb="irf" \
+           irf="cta_dummy_irf" \
+           outcube="bkgcube2.fits" \
+           outmodel="bkgcube2.xml"
+$ECHO -n "."
+if [ -s "bkgcube2.fits" ]
+then
+$ECHO -n "."
+else
+$ECHO " bkgcube2.fits file is not found"
+exit 1
+fi
+if [ -s "bkgcube2.fits" ]
+then
+$ECHO -n "."
+else
+$ECHO " bkgcube2.xml file is not found"
+exit 1
+fi
+$ECHO " ok"
+
 #
@@ -201,11 +342,12 @@ $ECHO &quot; ok&quot;
 $ECHO -n "Test ctmodel: "
 #
 # Run 1
-$ctmodel inobs="cntmap.fits" \
-         incube="cntmap.fits" \
+$ctmodel inobs="cntmap1.fits" \
+         incube="cntmap1.fits" \
          outcube="modmap1.fits" \
          expcube="NONE" \
          psfcube="NONE" \
+         bkgcube="NONE" \
          caldb="irf" \
          irf="cta_dummy_irf" \
          inmodel="data/crab.xml"
@@ -217,13 +359,35 @@ else
   $ECHO " modmap1.fits file is not found"
   exit 1
 fi
+
 #
 # Run 2
+$ctmodel inobs="cntmap2.fits" \
+         incube="cntmap2.fits" \
+         outcube="modmap2.fits" \
+         expcube="expcube2.fits" \
+         psfcube="psfcube2.fits" \
+         bkgcube="bkgcube2.fits" \
+         caldb="irf" \
+         irf="cta_dummy_irf" \
+         inmodel="bkgcube2.xml"
+$ECHO -n "."
+if [ -s "modmap2.fits" ]
+then
+$ECHO -n "."
+else
+$ECHO " modmap2.fits file is not found"
+exit 1
+fi
+
+#
+# Run 3
 $ctmodel inobs="NONE" \
          incube="NONE" \
          outcube="modmap3.fits" \
          expcube="NONE" \
          psfcube="NONE" \
+         bkgcube="NONE" \
          caldb="irf" \
          irf="cta_dummy_irf" \
          inmodel="data/crab.xml" \
@@ -250,24 +414,6 @@ else
   $ECHO " modmap3.fits file is not found"
   exit 1
 fi
-#
-# Run 3
-$ctmodel inobs="cntmap2.fits" \
-         incube="cntmap2.fits" \
-         outcube="modmap2.fits" \
-         expcube="NONE" \
-         psfcube="NONE" \
-         caldb="irf" \
-         irf="cta_dummy_irf" \
-         inmodel="data/model_background_cube.xml"
-$ECHO -n "."
-if [ -s "modmap2.fits" ]
-then
-  $ECHO -n "."
-else
-  $ECHO " modmap2.fits file is not found"
-  exit 1
-fi
 $ECHO " ok"
@@ -301,11 +447,12 @@ $ECHO &quot; ok&quot;
 $ECHO -n "Test ctlike: "
 #
 # Run 1
-$ctlike inobs="cntmap.fits" \
+$ctlike inobs="cntmap1.fits" \
         inmodel="data/crab.xml" \
         outmodel="results_binned.xml" \
         expcube="NONE" \
         psfcube="NONE" \
+        bkgcube="NONE" \
         caldb="irf" \
         irf="cta_dummy_irf"
 $ECHO -n "."
@@ -331,36 +478,23 @@ else
   $ECHO " results_unbinned.xml file is not found"
   exit 1
 fi
+
 #
 # Run 3
-$ctlike inobs="events2.fits" \
-        inmodel="data/model_background_cube.xml" \
-        outmodel="results_unbinned_background_cube.xml" \
-        caldb="irf" \
-        irf="cta_dummy_irf"
-$ECHO -n "."
-if [ -s "results_unbinned_background_cube.xml" ]
-then
-  $ECHO -n "."
-else
-  $ECHO " results_unbinned_background_cube.xml file is not found"
-  exit 1
-fi
-#
-# Run 4
 $ctlike inobs="cntmap2.fits" \
-        inmodel="data/model_background_cube.xml" \
-        outmodel="results_binned_background_cube.xml" \
-        expcube="NONE" \
-        psfcube="NONE" \
+        inmodel="bkgcube2.xml" \
+        outmodel="results_binned_cube_background.xml" \
+        expcube="expcube2.fits" \
+        psfcube="psfcube2.fits" \
+        bkgcube="bkgcube2.fits" \
         caldb="irf" \
         irf="cta_dummy_irf"
 $ECHO -n "."
-if [ -s "results_binned_background_cube.xml" ]
+if [ -s "results_binned_cube_background.xml" ]
 then
   $ECHO -n "."
 else
-  $ECHO " results_binned_background_cube.xml file is not found"
+  $ECHO " results_binned_cube_background.xml file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -394,167 +528,6 @@ fi
 $ECHO " ok"
-#
-# Test ctexpcube
-# ==============
-$ECHO -n "Test ctexpcube: "
-#
-# Run 1
-$ctexpcube inobs="data/crab_events.fits" \
-           incube="NONE" \
-           outcube="expcube1.fits" \
-           caldb="irf" \
-           irf="cta_dummy_irf" \
-           emin=0.1 \
-           emax=100.0 \
-           enumbins=20 \
-           nxpix=200 \
-           nypix=200 \
-           binsz=0.02 \
-           coordsys="CEL" \
-           xref=83.63 \
-           yref=22.01 \
-           proj="CAR"
-$ECHO -n "."
-if [ -s "expcube1.fits" ]
-then
-  $ECHO -n "."
-else
-  $ECHO " expcube1.fits file is not found"
-  exit 1
-fi
-#
-# Run 2
-$ctexpcube inobs="data/crab_events.fits" \
-           incube="data/crab_cntmap.fits" \
-           outcube="expcube2.fits" \
-           caldb="irf" \
-           irf="cta_dummy_irf"
-$ECHO -n "."
-if [ -s "expcube2.fits" ]
-then
-  $ECHO -n "."
-else
-  $ECHO " expcube2.fits file is not found"
-  exit 1
-fi
-$ECHO " ok"
-
-
-#
-# Test ctpsfcube
-# ==============
-$ECHO -n "Test ctpsfcube: "
-#
-# Run 1
-$ctpsfcube inobs="data/crab_events.fits" \
-           incube="NONE" \
-           outcube="psfcube1.fits" \
-           caldb="irf" \
-           irf="cta_dummy_irf" \
-           emin=0.1 \
-           emax=100.0 \
-           enumbins=20 \
-           nxpix=10 \
-           nypix=10 \
-           binsz=0.4 \
-           coordsys="CEL" \
-           xref=83.63 \
-           yref=22.01 \
-           proj="CAR" \
-           amax=0.3 \
-           anumbins=10
-$ECHO -n "."
-if [ -s "psfcube1.fits" ]
-then
-  $ECHO -n "."
-else
-  $ECHO " psfcube1.fits file is not found"
-  exit 1
-fi
-#
-# Run 2
-$ctpsfcube inobs="data/crab_events.fits" \
-           incube="data/crab_cntmap.fits" \
-           outcube="psfcube2.fits" \
-           caldb="irf" \
-           irf="cta_dummy_irf" \
-           amax=0.3 \
-           anumbins=2
-$ECHO -n "."
-if [ -s "psfcube2.fits" ]
-then
-  $ECHO -n "."
-else
-  $ECHO " psfcube2.fits file is not found"
-  exit 1
-fi
-$ECHO " ok"
-
-
-#
-# Test ctbkgcube
-# ==============
-$ECHO -n "Test ctbkgcube: "
-#
-# Run 1
-$ctbkgcube inobs="data/crab_events.fits" \
-           inmodel="data/crab.xml" \
-           incube="NONE" \
-           outcube="bkgcube1.fits" \
-           outmodel="bkgcube1.xml" \
-           caldb="irf" \
-           irf="cta_dummy_irf" \
-           emin=0.1 \
-           emax=100.0 \
-           enumbins=20 \
-           nxpix=10 \
-           nypix=10 \
-           binsz=0.4 \
-           coordsys="CEL" \
-           xref=83.63 \
-           yref=22.01 \
-           proj="CAR"
-$ECHO -n "."
-if [ -s "bkgcube1.fits" ]
-then
-  $ECHO -n "."
-else
-  $ECHO " bkgcube1.fits file is not found"
-  exit 1
-fi
-if [ -s "bkgcube1.xml" ]
-then
-  $ECHO -n "."
-else
-  $ECHO " bkgcube1.xml file is not found"
-  exit 1
-fi
-#
-# Run 2
-$ctbkgcube inobs="data/crab_events.fits" \
-           inmodel="data/crab.xml" \
-           incube="data/crab_cntmap.fits" \
-           caldb="irf" \
-           irf="cta_dummy_irf" \
-           outcube="bkgcube2.fits" \
-           outmodel="bkgcube2.xml"
-$ECHO -n "."
-if [ -s "bkgcube2.fits" ]
-then
-  $ECHO -n "."
-else
-  $ECHO " bkgcube2.fits file is not found"
-  exit 1
-fi
-if [ -s "bkgcube2.fits" ]
-then
-  $ECHO -n "."
-else
-  $ECHO " bkgcube2.xml file is not found"
-  exit 1
-fi
-$ECHO " ok"
 #
@@ -608,10 +581,10 @@ $ECHO &quot; ok&quot;
 # ===============
 $ECHO -n "Test ctulimit: "
 $ctulimit inobs="data/crab_events.fits.gz" \
-    inmodel="data/crab.xml" \
-    srcname="Crab" \
-    caldb="irf" \
-    irf="cta_dummy_irf"
+          inmodel="data/crab.xml" \
+          srcname="Crab" \
+          caldb="irf" \
+          irf="cta_dummy_irf"
 $ECHO -n "."
 if [ -s "ctulimit.log" ]
 then