Added support of 'edispcube' parameter to all ctools

Michael Mayer
1 parent b6c5a749
Showing 23 changed files with 144 additions and 122 deletions
cscripts/cslightcrv.py
cscripts/cspull.py
cscripts/csresmap.py
cscripts/csspec.py
cscripts/cstsdist.py
src/ctbutterfly/ctbutterfly.par
src/ctedispcube/ctedispcube.cpp
src/ctedispcube/ctedispcube.hpp
src/ctedispcube/ctedispcube.par
src/cterror/cterror.par
src/ctexpcube/ctexpcube.cpp
src/ctexpcube/ctexpcube.hpp
src/ctexpcube/ctexpcube.par
src/ctlike/ctlike.par
src/ctmodel/ctmodel.par
src/ctpsfcube/ctpsfcube.cpp
src/ctpsfcube/ctpsfcube.hpp
src/ctpsfcube/ctpsfcube.par
src/cttsmap/cttsmap.par
src/ctulimit/ctulimit.par
@@ -104,7 +104,9 @@ class cslightcrv(ctools.cscript):
             pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name"))
             pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Input exposure cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","Input PSF cube file (only needed for stacked analysis)"))
+            pars.append(gammalib.GApplicationPar("edispcube","f","a","NONE","","","Input energy dispersion cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Input background cube file (only needed for stacked analysis)"))
+            pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
             pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database"))
             pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function"))
             pars.append(gammalib.GApplicationPar("outfile","f","a","lightcurve.fits","","","Output light curve file"))
@@ -190,6 +192,9 @@ class cslightcrv(ctools.cscript):
         # Get source name   
         self.m_srcname = self["srcname"].string()
+        # Get energy dispersion flag
+        self.m_edisp = self["edisp"].boolean()
+        
         # Get time boundaries             
         self.create_tbounds()
@@ -472,12 +477,35 @@ class cslightcrv(ctools.cscript):
                 # Set new binned observation
                 obs = bin.obs()
-
+                
+                # Get new models
+                models = bkgcube.models()
+                
                 # Set precomputed binned response
                 obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.bkgcube())
-                # Get new models
-                models = bkgcube.models()
+#                 # Check if we need to include energy dispersion
+#                 if self.m_edisp:
+#                     
+#                     # Create edisp cube
+#                     edispcube = ctools.ctedispcube(select.obs())
+#                     edispcube["incube"]   = "NONE"
+#                     edispcube["usepnt"]   = False
+#                     edispcube["ebinalg"]  = "LOG"
+#                     edispcube["xref"]     = self.m_xref
+#                     edispcube["yref"]     = self.m_yref
+#                     edispcube["binsz"]    = self.m_binsz
+#                     edispcube["nxpix"]    = self.m_nxpix
+#                     edispcube["nypix"]    = self.m_nypix
+#                     edispcube["enumbins"] = self.m_ebins
+#                     edispcube["emin"]     = self.m_emin
+#                     edispcube["emax"]     = self.m_emax    
+#                     edispcube["coordsys"] = self.m_coordsys
+#                     edispcube["proj"]     = self.m_proj               
+#                     edispcube.run()
+#                     
+#                     # Set precomputed binned response including energy dispersion
+#                     obs[0].response(expcube.expcube(), psfcube.psfcube(), edispcube.edispcube(), bkgcube.bkgcube())                    
                 # Fix background models if required
                 if self.m_fix_bkg:
@@ -495,6 +523,7 @@ class cslightcrv(ctools.cscript):
             # Likelihood
             like = ctools.ctlike(obs)
+            like["edisp"] = self.m_edisp
             like.run()
             # Skip bin if no event was present
@@ -52,6 +52,7 @@ class cspull(ctools.cscript):
         self.m_exposure    = None
         self.m_psfcube     = None
         self.m_bckcube     = None
+        self.m_edispcube   = None
         self.m_stackmodels = None
         self.m_coordsys    = "CEL"
         self.m_proj        = "TAN"
@@ -246,6 +247,30 @@ class cspull(ctools.cscript):
         psfcube["proj"]     = self.m_proj
         psfcube.run()
+        if self.m_edisp:
+            
+            # Get stacked Edisp
+            edispcube = ctools.ctpsfcube(self.obs)
+            edispcube["incube"]   = "NONE"
+            edispcube["usepnt"]   = True
+            edispcube["ebinalg"]  = "LOG"
+            edispcube["binsz"]    = self.m_binsz*10.0
+            edispcube["nxpix"]    = self.m_npix/10
+            edispcube["nypix"]    = self.m_npix/10
+            edispcube["enumbins"] = self.m_enumbins
+            edispcube["emin"]     = self["emin"].real()
+            edispcube["emax"]     = self["emax"].real()
+            edispcube["coordsys"] = self.m_coordsys
+            edispcube["proj"]     = self.m_proj
+            edispcube.run()   
+            
+            # Store result
+            self.m_edispcube = edispcube.edispcube().copy()
+            
+            # Logging
+            if self.logTerse():
+                self.log("Computed Edisp cube\n")
+
         # Notify Psf computation
         if self.logTerse():
             self.log("Computed Psf cube\n")
@@ -391,7 +416,10 @@ class cspull(ctools.cscript):
         # If stacked, add stacked responses and model
         if self.obs.size() > 1 and self.m_enumbins > 0:
-            obs[0].response(self.m_exposure, self.m_psfcube, self.m_bckcube)
+            if self.m_edisp:
+                obs[0].response(self.m_exposure, self.m_psfcube, self.m_edispcube, self.m_bckcube)
+            else:
+                obs[0].response(self.m_exposure, self.m_psfcube, self.m_bckcube)
             obs.models(self.m_stackmodels)
         # Determine number of events in simulation
@@ -108,6 +108,7 @@ class csresmap(ctools.cscript):
             pars.append(gammalib.GApplicationPar("modcube","f","a","NONE","","","Input model cube file (generated with ctmodel)"))            
             pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Input exposure cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","Input PSF cube file (only needed for stacked analysis)"))
+            pars.append(gammalib.GApplicationPar("edispcube","f","a","NONE","","","Input energy dispersion cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("bkgcube","f","a","NONE","","","Input background cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Input model XML file"))
             pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database"))
@@ -105,6 +105,7 @@ class csspec(ctools.cscript):
             pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name"))
             pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Input exposure cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","Input PSF cube file (only needed for stacked analysis)"))
+            pars.append(gammalib.GApplicationPar("edispcube","f","a","NONE","","","Input energy dispersion cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("bkgcube","f","a","NONE","","","Input background cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database"))
             pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function"))
@@ -99,6 +99,7 @@ class cstsdist(ctools.cscript):
             pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name"))
             pars.append(gammalib.GApplicationPar("expcube","s","a","NONE","","","Input exposure cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("psfcube","s","a","NONE","","","Input PSF cube file (only needed for stacked analysis)"))
+            pars.append(gammalib.GApplicationPar("edispcube","f","a","NONE","","","Input energy dispersion cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Input background cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database"))
             pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function"))
@@ -22,6 +22,7 @@ inmodel,    f, a, &quot;$CTOOLS/share/models/crab.xml&quot;,,, &quot;Input model XML file&quot;
 srcname,    s, a, "Crab",,,"Source of interest"
 expcube,    f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
 psfcube,    f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
+edispcube,  f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)"
 bkgcube,    f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
 caldb,      s, a, "prod2",,, "Calibration database"
 irf,        s, a, "South_0.5h",,, "Instrument response function"
@@ -216,18 +216,6 @@ void ctedispcube::run(void)
         log << std::endl;
     }
-    // Set energy dispersion flag for all CTA observations and save old
-    // values in save_edisp vector
-//    std::vector<bool> save_edisp;
-//    save_edisp.assign(m_obs.size(), false);
-//    for (int i = 0; i < m_obs.size(); ++i) {
-//        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-//        if (obs != NULL) {
-//            save_edisp[i] = obs->response()->apply_edisp();
-//            obs->response()->apply_edisp(m_apply_edisp);
-//        }
-//    }
-
     // Write observation(s) into logger
     if (logTerse()) {
         log << std::endl;
@@ -246,31 +234,23 @@ void ctedispcube::run(void)
         log.header1("Generate EDISP cube");
     }
-    // Fill EDISP
+    // Fill Edisp
     m_edispcube.fill(m_obs, &log);
-    // Restore energy dispersion flag for all CTA observations
-//    for (int i = 0; i < m_obs.size(); ++i) {
-//        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-//        if (obs != NULL) {
-//            obs->response()->apply_edisp(save_edisp[i]);
-//        }
-//    }
-
     // Return
     return;
 }
 /***********************************************************************//**
- * @brief Save EDISP cube
+ * @brief Save Edisp cube
  ***************************************************************************/
 void ctedispcube::save(void)
 {
     // Write header
     if (logTerse()) {
         log << std::endl;
-        log.header1("Save EDISP cube");
+        log.header1("Save Edisp cube");
     }
     // Get output filename
@@ -281,7 +261,7 @@ void ctedispcube::save(void)
         // Log filename
         if (logTerse()) {
-            log << "Save EDISP cube into file \""+m_outcube+"\".";
+            log << "Save Edisp cube into file \""+m_outcube+"\".";
             log << std::endl;
         }
@@ -308,7 +288,6 @@ void ctedispcube::init_members(void)
 {
     // Initialise members
     m_outcube.clear();
-//    m_apply_edisp = false;
     // Initialise protected members
     m_obs.clear();
@@ -328,7 +307,6 @@ void ctedispcube::copy_members(const ctedispcube&amp; app)
 {
     // Copy attributes
     m_outcube     = app.m_outcube;
-//    m_apply_edisp = app.m_apply_edisp;
     // Copy protected members
     m_obs        = app.m_obs;
@@ -396,14 +374,11 @@ void ctedispcube::get_parameters(void)
         // Load event cube from filename
         GCTAEventCube cube(incube);
-        // Define edisp cube
+        // Define Edisp cube
         m_edispcube = GCTACubeEdisp(cube, mmax, migrabins);
     } // endelse: cube loaded from file
-    // Read energy dispersion flag
-//    m_apply_edisp = (*this)["edisp"].boolean();
-
     // Read output filename (if needed)
     if (read_ahead()) {
         m_outcube = (*this)["outcube"].filename();
@@ -69,12 +69,11 @@ protected:
     void get_parameters(void);
     // User parameters
-    GFilename     m_outcube;     //!< Output EDISP cube file name
-    bool          m_apply_edisp; //!< Apply energy dispersion?
+    GFilename     m_outcube;     //!< Output Edisp cube file name
     // Protected members
     GObservations m_obs;         //!< Observation container
-    GCTACubeEdisp   m_edispcube;     //!< EDISP cube
+    GCTACubeEdisp   m_edispcube;     //!< Edisp cube
 };
@@ -21,8 +21,7 @@ inobs,    f, a, &quot;NONE&quot;,,, &quot;Input event list or observation definition XML file&quot;
 incube,   f, a, "NONE",,, "Input counts cube file to extract PSF cube definition"
 caldb,    s, a, "prod2",,, "Calibration database"
 irf,      s, a, "South_0.5h",,, "Instrument response function"
-edisp,    b, h, no,,, "Apply energy dispersion?"
-outcube,  f, a, "psfcube.fits",,, "Output PSF cube file"
+outcube,  f, a, "edispcube.fits",,, "Output PSF cube file"
 #
 # Energy binning parameters
@@ -58,4 +57,4 @@ chatter, i, h, 2,0,4, &quot;Chattiness of output&quot;
 clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
 debug,   b, h, no,,, "Debugging mode activated"
 mode,    s, h, "ql",,, "Mode of automatic parameters"
-logfile, f, h, "ctpsfcube.log",,, "Log filename"
+logfile, f, h, "ctedispcube.log",,, "Log filename"
@@ -23,6 +23,7 @@ inmodel,  f, a, &quot;$CTOOLS/share/models/crab.xml&quot;,,, &quot;Input model XML file&quot;
 srcname,  s, a, "Crab",,, "Source of interest"
 expcube,  f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
 psfcube,  f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
+edispcube,f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)"
 bkgcube,  f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
 caldb,    s, a, "prod2",,, "Calibration database"
 irf,      s, a, "South_0.5h",,, "Instrument response function"
@@ -217,18 +217,6 @@ void ctexpcube::run(void)
         log << std::endl;
     }
-    // Set energy dispersion flag for all CTA observations and save old
-    // values in save_edisp vector
-    std::vector<bool> save_edisp;
-    save_edisp.assign(m_obs.size(), false);
-    for (int i = 0; i < m_obs.size(); ++i) {
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-        if (obs != NULL) {
-            save_edisp[i] = obs->response()->apply_edisp();
-            obs->response()->apply_edisp(m_apply_edisp);
-        }
-    }
-
     // Write observation(s) into logger
     if (logTerse()) {
         log << std::endl;
@@ -250,14 +238,6 @@ void ctexpcube::run(void)
     // Fill exposure
     m_expcube.fill(m_obs, &log);
-    // Restore energy dispersion flag for all CTA observations
-    for (int i = 0; i < m_obs.size(); ++i) {
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-        if (obs != NULL) {
-            obs->response()->apply_edisp(save_edisp[i]);
-        }
-    }
-
     // Optionally publish exposure cube
     if (m_publish) {
         publish();
@@ -348,7 +328,6 @@ void ctexpcube::init_members(void)
 {
     // Initialise members
     m_outcube.clear();
-    m_apply_edisp = false;
     m_publish     = false;
     // Initialise protected members
@@ -369,7 +348,6 @@ void ctexpcube::copy_members(const ctexpcube&amp; app)
 {
     // Copy attributes
     m_outcube     = app.m_outcube;
-    m_apply_edisp = app.m_apply_edisp;
     m_publish     = app.m_publish;
     // Copy protected members
@@ -414,9 +392,6 @@ void ctexpcube::get_parameters(void)
     } // endif: there was no observation in the container
-    // Read energy dispersion flag
-    m_apply_edisp = (*this)["edisp"].boolean();
-
     // Get the incube filename
     std::string incube = (*this)["incube"].filename();
@@ -71,7 +71,6 @@ protected:
     // User parameters
     GFilename m_outcube;     //!< Output exposure cube file
-    bool      m_apply_edisp; //!< Apply energy dispersion?
     bool      m_publish;     //!< Publish exposure cube?
     // Protected members
@@ -21,7 +21,6 @@ inobs,    f, a, &quot;NONE&quot;,,, &quot;Input event list or observation definition XML file&quot;
 incube,   f, a, "NONE",,, "Input counts cube file to extract exposure cube definition"
 caldb,    s, a, "prod2",,, "Calibration database"
 irf,      s, a, "South_0.5h",,, "Instrument response function"
-edisp,    b, h, no,,, "Apply energy dispersion?"
 outcube,  f, a, "expcube.fits",,, "Output exposure cube file"
 #
@@ -21,6 +21,7 @@ inobs,           f, a, &quot;events.fits&quot;,,, &quot;Input event list, counts cube or observ
 inmodel,         f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file"
 expcube,         f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
 psfcube,         f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
+edispcube,       f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)"
 bkgcube,         f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
 caldb,           s, a, "prod2",,, "Calibration database"
 irf,             s, a, "South_0.5h",,, "Instrument response function"
@@ -22,6 +22,7 @@ incube,   f, a, &quot;NONE&quot;,,, &quot;Input counts cube file to extract model cube definiti
 inmodel,  f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file"
 expcube,  f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
 psfcube,  f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
+edispcube,f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)"
 bkgcube,  f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
 caldb,    s, a, "prod2",,, "Calibration database"
 irf,      s, a, "South_0.5h",,, "Instrument response function"
@@ -216,18 +216,6 @@ void ctpsfcube::run(void)
         log << std::endl;
     }
-    // Set energy dispersion flag for all CTA observations and save old
-    // values in save_edisp vector
-    std::vector<bool> save_edisp;
-    save_edisp.assign(m_obs.size(), false);
-    for (int i = 0; i < m_obs.size(); ++i) {
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-        if (obs != NULL) {
-            save_edisp[i] = obs->response()->apply_edisp();
-            obs->response()->apply_edisp(m_apply_edisp);
-        }
-    }
-
     // Write observation(s) into logger
     if (logTerse()) {
         log << std::endl;
@@ -249,14 +237,6 @@ void ctpsfcube::run(void)
     // Fill PSF 
     m_psfcube.fill(m_obs, &log);
-    // Restore energy dispersion flag for all CTA observations
-    for (int i = 0; i < m_obs.size(); ++i) {
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-        if (obs != NULL) {
-            obs->response()->apply_edisp(save_edisp[i]);
-        }
-    }
-
     // Return
     return;
 }
@@ -308,7 +288,6 @@ void ctpsfcube::init_members(void)
 {
     // Initialise members
     m_outcube.clear();
-    m_apply_edisp = false;
     // Initialise protected members
     m_obs.clear();
@@ -328,7 +307,6 @@ void ctpsfcube::copy_members(const ctpsfcube&amp; app)
 {
     // Copy attributes
     m_outcube     = app.m_outcube;
-    m_apply_edisp = app.m_apply_edisp;
     // Copy protected members
     m_obs        = app.m_obs;
@@ -401,9 +379,6 @@ void ctpsfcube::get_parameters(void)
     } // endelse: cube loaded from file
-    // Read energy dispersion flag
-    m_apply_edisp = (*this)["edisp"].boolean();
-
     // Read output filename (if needed)
     if (read_ahead()) {
         m_outcube = (*this)["outcube"].filename();
@@ -70,7 +70,6 @@ protected:
     // User parameters
     GFilename     m_outcube;     //!< Output PSF cube file name
-    bool          m_apply_edisp; //!< Apply energy dispersion?
     // Protected members
     GObservations m_obs;         //!< Observation container
@@ -21,7 +21,6 @@ inobs,    f, a, &quot;NONE&quot;,,, &quot;Input event list or observation definition XML file&quot;
 incube,   f, a, "NONE",,, "Input counts cube file to extract PSF cube definition"
 caldb,    s, a, "prod2",,, "Calibration database"
 irf,      s, a, "South_0.5h",,, "Instrument response function"
-edisp,    b, h, no,,, "Apply energy dispersion?"
 outcube,  f, a, "psfcube.fits",,, "Output PSF cube file"
 #
@@ -17,35 +17,36 @@
 #
 # General parameters
 #===================
-inobs,   f, a, "events.fits",,, "Input event list, counts cube or observation definition XML file"
-inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file"
-srcname, s, a, "Crab",,, "Test source name"
-expcube, f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
-psfcube, f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
-bkgcube, f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
-caldb,   s, a, "prod2",,, "Calibration database"
-irf,     s, a, "South_0.5h",,, "Instrument response function"
-edisp,    b, h, no,,, "Apply energy dispersion?"
-outmap,  f, a, "tsmap.fits",,, "Output Test Statistic map file"
+inobs,     f, a, "events.fits",,, "Input event list, counts cube or observation definition XML file"
+inmodel,   f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file"
+srcname,   s, a, "Crab",,, "Test source name"
+expcube,   f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
+psfcube,   f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
+edispcube, f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)"
+bkgcube,   f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
+caldb,     s, a, "prod2",,, "Calibration database"
+irf,       s, a, "South_0.5h",,, "Instrument response function"
+edisp,     b, h, no,,, "Apply energy dispersion?"
+outmap,    f, a, "tsmap.fits",,, "Output Test Statistic map file"
 #
 # Spatial binning parameters
 #===========================
-usepnt,   b, h, no,,, "Use pointing instead of xref/yref parameters?"
-nxpix,    i, a, 200,,, "Size of the X axis in pixels"
-nypix,    i, a, 200,,, "Size of the Y axis in pixels"
-binsz,    r, a, 0.02,,, "Image scale (in degrees/pixel)"
-coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
-xref,     r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
-yref,     r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
-proj,     s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
+usepnt,    b, h, no,,, "Use pointing instead of xref/yref parameters?"
+nxpix,     i, a, 200,,, "Size of the X axis in pixels"
+nypix,     i, a, 200,,, "Size of the Y axis in pixels"
+binsz,     r, a, 0.02,,, "Image scale (in degrees/pixel)"
+coordsys,  s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+xref,      r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
+yref,      r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
+proj,      s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
 #
 # Parameters for splitting and speed purpose
 #===========================
-binmin,   i, h, -1,,, "First bin to compute"
-binmax,   i, h, -1,,, "Last bin to compute"
-logL0,    r, h, 0.0,,, "LogLikelihood value of null hypothesis"
+binmin,    i, h, -1,,, "First bin to compute"
+binmax,    i, h, -1,,, "Last bin to compute"
+logL0,     r, h, 0.0,,, "LogLikelihood value of null hypothesis"
 #
 # Standard parameters
@@ -22,6 +22,7 @@ inmodel,  f, a, &quot;$CTOOLS/share/models/crab.xml&quot;,,, &quot;Input model XML file&quot;
 srcname,  s, a, "Crab",,, "Source of interest"
 expcube,  f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
 psfcube,  f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
+edispcube,f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)"
 bkgcube,  f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
 caldb,    s, a, "prod2",,, "Calibration database"
 irf,      s, a, "South_0.5h",,, "Instrument response function"
@@ -862,11 +862,12 @@ 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") && has_par("bkgcube")) {
+        if (has_par("expcube") && has_par("psfcube") && has_par("bkgcube") && has_par("edispcube")) {
             // Get filenames
             std::string expcube = (*this)["expcube"].filename();
             std::string psfcube = (*this)["psfcube"].filename();
+            std::string edispcube = (*this)["edispcube"].filename();
             std::string bkgcube = (*this)["bkgcube"].filename();
             // Extract response information if available
@@ -880,8 +881,15 @@ void ctool::set_obs_response(GCTAObservation* obs)
                 GCTACubePsf        psf(psfcube);
                 GCTACubeBackground background(bkgcube);
-                // Set reponse
-                obs->response(exposure, psf, background);
+                // Optionally load Edisp cube
+                if (edispcube != "NONE" && (gammalib::strip_whitespace(edispcube) != "")) {
+                	GCTACubeEdisp edisp(edispcube);
+                	obs->response(exposure, psf, edisp, background);
+                }
+                else {
+					// Set reponse
+					obs->response(exposure, psf, background);
+                }
                 // Signal response availability
                 has_response = true;
@@ -405,6 +405,7 @@ fi
   inobs="data/crab_cntmap_small.fits.gz" \
   expcube="data/crab_expcube.fits.gz" \
   psfcube="data/crab_psfcube.fits.gz" \
+  edispcube="data/crab_edispcube.fits.gz" \
   bkgcube="data/crab_bkgcube.fits.gz" \
   inmodel="data/crab_bkgcube.xml" \
   ebinalg="LOG" \
@@ -293,8 +293,8 @@ $ctedispcube inobs=&quot;data/crab_events.fits&quot; \
              xref=83.63 \
              yref=22.01 \
              proj="CAR" \
-             mmax=10.0 \
-             migrabins=100
+             mmax=3.0 \
+             migrabins=10
 $ECHO -n "."
 if [ -s "edispcube1.fits" ]
 then
@@ -310,8 +310,8 @@ $ctedispcube inobs=&quot;events.fits&quot; \
              outcube="edispcube2.fits" \
              caldb="irf" \
              irf="cta_dummy_irf" \
-             mmax=10.0 \
-             migrabins=5
+             mmax=3.0 \
+             migrabins=3
 $ECHO -n "."
 if [ -s "edispcube2.fits" ]
 then
@@ -401,6 +401,7 @@ $ctmodel inobs=&quot;cntmap1.fits&quot; \
          expcube="NONE" \
          psfcube="NONE" \
          bkgcube="NONE" \
+         edispcube="NONE" \
          caldb="irf" \
          irf="cta_dummy_irf" \
          inmodel="data/crab.xml"
@@ -420,6 +421,7 @@ $ctmodel inobs=&quot;cntmap2.fits&quot; \
          outcube="modmap2.fits" \
          expcube="expcube2.fits" \
          psfcube="psfcube2.fits" \
+         edispcube="edispcube2.fits" \
          bkgcube="bkgcube2.fits" \
          caldb="irf" \
          irf="cta_dummy_irf" \
@@ -440,6 +442,7 @@ $ctmodel inobs=&quot;NONE&quot; \
          outcube="modmap3.fits" \
          expcube="NONE" \
          psfcube="NONE" \
+         edispcube="NONE" \
          bkgcube="NONE" \
          caldb="irf" \
          irf="cta_dummy_irf" \
@@ -505,6 +508,7 @@ $ctlike inobs=&quot;cntmap1.fits&quot; \
         outmodel="results_binned.xml" \
         expcube="NONE" \
         psfcube="NONE" \
+        edispcube="NONE" \
         bkgcube="NONE" \
         caldb="irf" \
         irf="cta_dummy_irf"
@@ -539,6 +543,7 @@ $ctlike inobs=&quot;cntmap2.fits&quot; \
         outmodel="results_binned_cube_background.xml" \
         expcube="expcube2.fits" \
         psfcube="psfcube2.fits" \
+        edispcube="NONE" \
         bkgcube="bkgcube2.fits" \
         caldb="irf" \
         irf="cta_dummy_irf"
@@ -550,9 +555,32 @@ else
   $ECHO " results_binned_cube_background.xml file is not found"
   exit 1
 fi
+
+
+#
+# Run 4
+$ctlike inobs="cntmap2.fits" \
+        inmodel="bkgcube2.xml" \
+        outmodel="results_binned_cube_background_edisp.xml" \
+        expcube="expcube2.fits" \
+        psfcube="psfcube2.fits" \
+        edispcube="edispcube2.fits" \
+        bkgcube="bkgcube2.fits" \
+        edisp="yes"
+        caldb="irf" \
+        irf="cta_dummy_irf"
+$ECHO -n "."
+if [ -s "results_binned_cube_background_edisp.xml" ]
+then
+    $ECHO -n "."
+else
+    $ECHO " results_binned_cube_background_edisp.xml file is not found"
+    exit 1
+fi
 $ECHO " ok"
+
 #
 # Test cttsmap
 # ============