added support of edisp parameter to ctools and cscripts

mimayer
1 parent 8caa108d
Showing 16 changed files with 160 additions and 6 deletions
cscripts/cshessobs.py
cscripts/csresmap.py
cscripts/cssens.py
cscripts/csspec.py
src/ctbutterfly/ctbutterfly.cpp
src/ctbutterfly/ctbutterfly.hpp
src/ctbutterfly/ctbutterfly.par
src/cterror/cterror.cpp
src/cterror/cterror.hpp
src/cterror/cterror.par
src/cttsmap/cttsmap.cpp
src/cttsmap/cttsmap.hpp
src/cttsmap/cttsmap.par
src/ctulimit/ctulimit.cpp
src/ctulimit/ctulimit.hpp
src/ctulimit/ctulimit.par
@@ -45,6 +45,7 @@ class cshessobs(ctools.cscript):
         self.name    = "cshessobs"
         self.version = "0.1.0"
         self.datapath = ""
+        self.m_erange = gammalib.GEbounds()
         # Check on existence of HESSFITS environment variable
         try:
@@ -320,6 +321,12 @@ class cshessobs(ctools.cscript):
         # Return model
         return bck
+    
+    def erange(self):
+        """
+        Returns runlist energy range
+        """
+        return self.m_erange
     def execute(self):
         """
@@ -397,6 +404,10 @@ class cshessobs(ctools.cscript):
             self.log("\n")
             self.log.header1("Looping over runs")
+        # Initialise energy range values for logging
+        runlist_emin = 100.0
+        runlist_emax = 0.0
+
         for runnr in runlist:
             # Unify run number string length
@@ -411,7 +422,7 @@ class cshessobs(ctools.cscript):
             eventfile = runfolder + "hess_events_"+run+".fits.gz"
             aefffile  = runfolder + "hess_aeff_2d_"+run+".fits.gz"
             psffile   = runfolder + "hess_psf_king_"+run+".fits.gz"
-            edispfile = ""#runfolder + "hess_edisp_2d_"+run+".fits.gz"
+            edispfile = runfolder + "hess_edisp_2d_"+run+".fits.gz"
             bgfile    = runfolder + "hess_bkg_offruns_"+run+".fits.gz"
             # Check for existence of files
@@ -429,6 +440,9 @@ class cshessobs(ctools.cscript):
             elif not os.path.isfile(psffile):
                 msg = "Run "+str(int(run))+" has no PSF - Run is skipped"
                 skip = True
+            elif not os.path.isfile(edispfile):
+                msg = "Run "+str(int(run))+" has no energy dispersion - Run is skipped"
+                skip = True
             elif not os.path.isfile(bgfile):
                 msg = "Run "+str(int(run))+" has no background - Run is skipped"
                 skip = True            
@@ -458,6 +472,16 @@ class cshessobs(ctools.cscript):
             # Close FITS file
             fits.close()
+            
+            # Open effective area file to look for threshold
+            aeff_fits = gammalib.GFits(aefffile)
+            run_emin = aeff_fits["EFFECTIVE AREA"].real("LO_THRES")
+            run_emax = aeff_fits["EFFECTIVE AREA"].real("HI_THRES")
+            if run_emin < runlist_emin:
+                runlist_emin = run_emin
+            if run_emax > runlist_emax:
+                runlist_emax = run_emax
+            aeff_fits.close()          
             # Append observation to XML and set attributes
             obs = lib.append("observation");
@@ -494,6 +518,15 @@ class cshessobs(ctools.cscript):
             # Append instrumental background model
             self.models.append(self.hess_inst_background(run,aefffile,trgrate))
+        
+        # Store energy range as member
+        self.m_erange = gammalib.GEbounds(gammalib.GEnergy(runlist_emin,"TeV"),gammalib.GEnergy(runlist_emax,"TeV"))
+           
+        # Write energy range into log file    
+        if self.logTerse():
+            self.log("\n")
+            self.log("Eenrgy range of obervation list: ")
+            self.log(("% 3.2f" % runlist_emin)+" - "+("% 3.2f" % runlist_emax)+" TeV")
         # Append models provided by 'inmodels' of necessary
         if not self.inmodels == None:
@@ -110,6 +110,7 @@ class csresmap(ctools.cscript):
             pars.append(gammalib.GApplicationPar("bkgcube","f","a","NONE","","","Background cube file (only needed for stacked analysis)"))
             pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))
             pars.append(gammalib.GApplicationPar("outmap","f","a","resmap.fits","","","Output residual map"))
+            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_50h","","","Instrument response function"))
             pars.append(gammalib.GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm"))
@@ -141,7 +142,7 @@ class csresmap(ctools.cscript):
         self.m_skip_binning = False
         # First check if the inobs parameter is a counts cube
-        if self["inobs"].filename() != "NONE":
+        if self.obs.size() == 0 and self["inobs"].filename() != "NONE":
             if gammalib.is_fits(self["inobs"].filename()):
                 cta = gammalib.GCTAObservation()
                 cta.load(self["inobs"].filename())
@@ -189,6 +190,7 @@ class csresmap(ctools.cscript):
                 self.m_nxpix     = self["nxpix"].integer()
                 self.m_nypix     = self["nypix"].integer()
                 self.m_binsz     = self["binsz"].real()
+                self.m_edisp     = self["edisp"].boolean()
         # Read necessary parameters
         self.m_outfile   = self["outmap"].filename()    
@@ -301,6 +303,7 @@ class csresmap(ctools.cscript):
             model["chatter"].integer(self.m_chatter)
             model["clobber"].boolean(self.m_clobber)
             model["debug"].boolean(self.m_debug)
+            model["edisp"].boolean(self.m_edisp)
             model.run()
             # Get model map into GSkyMap object
@@ -107,6 +107,7 @@ class cssens(ctools.cscript):
             pars.append(gammalib.GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name"))
             pars.append(gammalib.GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
             pars.append(gammalib.GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)"))
+            pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
             pars.append(gammalib.GApplicationPar("emin","r","a","0.020","","","Lower energy limit (TeV)"))
             pars.append(gammalib.GApplicationPar("emax","r","a","200.0","","","Upper energy limit (TeV)"))
             pars.append(gammalib.GApplicationPar("bins","i","a","21","","","Number of energy bins for differential sensitivity computation"))
@@ -156,6 +157,7 @@ class cssens(ctools.cscript):
             self.m_binsz = 0.0
         # Read remaining parameters
+        self.m_edisp    = self["edisp"].boolean()
         self.m_ts_thres = self["sigma"].real()*self["sigma"].real()
         self.m_max_iter = self["max_iter"].integer()
         self.m_num_avg  = self["num_avg"].integer()
@@ -466,7 +468,7 @@ class cssens(ctools.cscript):
             # Simulate events
             sim = obsutils.sim(obs, nbins=self.m_enumbins, seed=iter,
                                binsz=self.m_binsz, npix=self.m_npix,
-                               log=self.m_log, debug=self.m_debug)
+                               log=self.m_log, debug=self.m_debug, edisp=self.m_edisp)
             # Determine number of events in simulation
             nevents = 0.0
@@ -482,7 +484,7 @@ class cssens(ctools.cscript):
             # Fit background only
             sim.models(bkg_model)
-            like       = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
+            like       = obsutils.fit(sim, log=self.m_log, debug=self.m_debug, edisp=self.m_edisp)
             result_bgm = like.obs().models().copy()
             LogL_bgm   = like.opt().value()
             npred_bgm  = like.obs().npred()
@@ -520,7 +522,7 @@ class cssens(ctools.cscript):
             # Fit background and test source
             sim.models(src_model)
-            like       = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
+            like       = obsutils.fit(sim, log=self.m_log, debug=self.m_debug, edisp=self.m_edisp)
             result_all = like.obs().models().copy()
             LogL_all   = like.opt().value()
             npred_all  = like.obs().npred()
@@ -107,6 +107,7 @@ class csspec(ctools.cscript):
             pars.append(gammalib.GApplicationPar("psfcube","f","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","prod2","","","Calibration database"))
+            pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
             pars.append(gammalib.GApplicationPar("irf","s","a","South_50h","","","Instrument response function"))
             pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit for spectral points(TeV)"))
             pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit for spectral points(TeV)"))
@@ -153,6 +154,10 @@ class csspec(ctools.cscript):
         # Get binning flag
         self.m_binned = self["binned"].boolean()
+        
+        # Get edisp flag
+        self.m_edisp = self["edisp"].boolean()
+        
         if self.m_binned:
             self.m_xref     = self["xref"].real()
             self.m_yref     = self["yref"].real()
@@ -368,6 +373,7 @@ class csspec(ctools.cscript):
                 expcube["incube"]   = "NONE"
                 expcube["usepnt"]   = False
                 expcube["ebinalg"]  = "LOG"
+                expcube["edisp"]    = self.m_edisp
                 expcube["xref"]     = self.m_xref
                 expcube["yref"]     = self.m_yref
                 expcube["binsz"]    = self.m_binsz
@@ -389,6 +395,7 @@ class csspec(ctools.cscript):
                 psfcube["incube"]   = "NONE"
                 psfcube["usepnt"]   = False
                 psfcube["ebinalg"]  = "LOG"
+                psfcube["edisp"]    = self.m_edisp
                 psfcube["xref"]     = self.m_xref
                 psfcube["yref"]     = self.m_yref
                 psfcube["binsz"]    = self.m_binsz
@@ -447,6 +454,7 @@ class csspec(ctools.cscript):
             # Likelihood
             like = ctools.ctlike(obs)
+            like["edisp"].boolean(self.m_edisp)
             like.run()
             # Skip bin if no event was present
@@ -216,6 +216,18 @@ void ctbutterfly::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;
@@ -347,6 +359,14 @@ void ctbutterfly::run(void)
     } //endfor: loop over energy bin
+    // 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;
 }
@@ -402,6 +422,7 @@ void ctbutterfly::init_members(void)
     m_srcname.clear();
     m_outfile.clear();
     m_ebounds.clear();
+    m_apply_edisp = false;
     // Initialise protected members
     m_obs.clear();
@@ -426,6 +447,7 @@ void ctbutterfly::copy_members(const ctbutterfly&amp; app)
     m_srcname = app.m_srcname;
     m_outfile = app.m_outfile;
     m_ebounds = app.m_ebounds;
+    m_apply_edisp = app.m_apply_edisp;
     // Copy protected members
     m_obs        = app.m_obs;
@@ -511,6 +533,9 @@ void ctbutterfly::get_parameters(void)
         // m_covariance.load(matrixfilename);
     }
+    // Read energy dispersion flag
+    m_apply_edisp = (*this)["edisp"].boolean();
+
     // Optionally read ahead parameters so that they get correctly
     // dumped into the log file
     if (read_ahead()) {
@@ -83,6 +83,7 @@ protected:
     std::string m_srcname;    //!< Name of source to compute butterfly
     std::string m_outfile;    //!< Output ascii file
     GEbounds    m_ebounds;    //!< Energy binning definition
+    bool          m_apply_edisp; //!< Apply energy dispersion?
     // Protected members
     GObservations       m_obs;        //!< Observation container
@@ -27,6 +27,7 @@ caldb,    s, a, &quot;prod2&quot;,,, &quot;Calibration database&quot;
 irf,      s, a, "South_50h",,, "Instrument response function"
 outfile,  f, a, "butterfly.txt",,, "Output ascii file"
 matrix,   f, h, "NONE",,, "Input covariance Matrix FITS file"
+edisp,    b, h, no,,, "Apply energy dispersion?"
 #
 # Energy binning parameters
@@ -214,6 +214,18 @@ void cterror::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;
@@ -355,6 +367,14 @@ void cterror::run(void)
     // Recover optimizer
     m_opt = best_opt;
+    // 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;
 }
@@ -410,6 +430,7 @@ void cterror::init_members(void)
     m_dlogL        = 0.0;
     m_best_logL    = 0.0;
     m_model_par    = NULL;
+    m_apply_edisp = false;
     // Set optimizer parameters
     m_opt.max_iter(m_max_iter);
@@ -432,6 +453,7 @@ void cterror::copy_members(const cterror&amp; app)
     m_confidence = app.m_confidence;
     m_tol        = app.m_tol;
     m_max_iter   = app.m_max_iter;
+    m_apply_edisp = app.m_apply_edisp;
     // Copy protected members
     m_obs       = app.m_obs;
@@ -500,6 +522,9 @@ void cterror::get_parameters(void)
         throw GException::invalid_value(G_GET_PARAMETERS, msg);
     }
+    // Read energy dispersion flag
+    m_apply_edisp = (*this)["edisp"].boolean();
+
     // Get confidence level and transform into log-likelihood difference
     m_confidence = (*this)["confidence"].real();
     double sigma = gammalib::erfinv(m_confidence) * gammalib::sqrt_two;
@@ -78,6 +78,7 @@ protected:
     double        m_tol;          //!< Tolerance for limit determination
     int           m_max_iter;     //!< Maximum number of iterations
     double        m_value;        //!< Parameter value 
+    bool          m_apply_edisp; //!< Apply energy dispersion?
     // Protected members
     GObservations m_obs;          //!< Observation container
+
 #############################################################################
 #                                                                           #
 #                       CTA Science Analysis Tools                          #
@@ -26,6 +27,7 @@ psfcube,  f, a, &quot;NONE&quot;,,, &quot;PSF cube file (only needed for stacked analysis)&quot;
 bkgcube,  f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
 caldb,    s, a, "prod2",,, "Calibration database"
 irf,      s, a, "South_50h",,, "Instrument response function"
+edisp,    b, h, no,,, "Apply energy dispersion?"
 #
 # Error calculation parameters
@@ -215,6 +215,18 @@ void cttsmap::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;
@@ -348,6 +360,14 @@ void cttsmap::run(void)
     // Bring models to initial state
     m_obs.models(models_orig);
+    // 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;
 }
@@ -411,6 +431,7 @@ void cttsmap::init_members(void)
     // Initialise members
     m_srcname.clear();
     m_outmap.clear();
+    m_apply_edisp = false;
     // Initialise protected members
     m_obs.clear();
@@ -438,6 +459,7 @@ void cttsmap::copy_members(const cttsmap&amp; app)
     // Copy attributes
     m_srcname = app.m_srcname;
     m_outmap  = app.m_outmap;
+    m_apply_edisp = app.m_apply_edisp;
     // Copy protected members
     m_binmin    = app.m_binmin;
@@ -542,6 +564,9 @@ void cttsmap::get_parameters(void)
     // Initialise maps from user parameters
     init_maps(map);
+    // Read energy dispersion flag
+    m_apply_edisp = (*this)["edisp"].boolean();
+
     // Get optional splitting parameters
     m_binmin = (*this)["binmin"].integer();
     m_binmax = (*this)["binmax"].integer();
@@ -86,6 +86,7 @@ protected:
     // User parameters
     std::string              m_srcname;    //!< Name of source which is moved around
     std::string              m_outmap;     //!< Output counts map or XML file
+    bool          m_apply_edisp; //!< Apply energy dispersion?
     // Parameters to control speed and job splitting
     int                      m_binmin;     //!< Map bin number from which computation should start
@@ -26,6 +26,7 @@ bkgcube, f, a, &quot;NONE&quot;,,, &quot;Background cube file (only needed for stacked analysis
 caldb,   s, a, "prod2",,, "Calibration database"
 irf,     s, a, "South_50h",,, "Instrument response function"
 outmap,  f, a, "tsmap.fits",,, "Output Test Statistic map"
+edisp,    b, h, no,,, "Apply energy dispersion?"
 #
 # Spatial binning parameters
@@ -217,6 +217,18 @@ void ctulimit::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;
@@ -343,6 +355,14 @@ void ctulimit::run(void)
     // Recover optimizer
     m_opt = best_opt;
+    // 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;
 }
@@ -405,6 +425,7 @@ void ctulimit::init_members(void)
     m_emax       = 0.0;
     m_tol        = 1.0e-6;
     m_max_iter   = 50;
+    m_apply_edisp = false;
     // Initialise protected members
     m_obs.clear();
@@ -443,7 +464,7 @@ void ctulimit::copy_members(const ctulimit&amp; app)
     m_emax       = app.m_emax;
     m_tol        = app.m_tol;
     m_max_iter   = app.m_max_iter;
-
+    m_apply_edisp = app.m_apply_edisp;
     // Copy protected members
     m_obs          = app.m_obs;
@@ -520,6 +541,9 @@ void ctulimit::get_parameters(void)
     // Get relevant model and parameter for upper limit computation
     get_model_parameter();
+    // Read energy dispersion flag
+    m_apply_edisp = (*this)["edisp"].boolean();
+
     // Get confidence level and transform into log-likelihood difference
     m_confidence = (*this)["confidence"].real();
     double sigma = gammalib::erfinv(m_confidence) * gammalib::sqrt_two;
@@ -93,6 +93,7 @@ protected:
     double        m_emax;         //!< Maximum energy for flux limits (TeV)
     double        m_tol;          //!< Tolerance for limit determination
     int           m_max_iter;     //!< Maximum number of iterations
+    bool          m_apply_edisp; //!< Apply energy dispersion?
     // Protected members
     GObservations m_obs;          //!< Observation container
@@ -25,6 +25,7 @@ psfcube,  f, a, &quot;NONE&quot;,,, &quot;PSF cube file (only needed for stacked analysis)&quot;
 bkgcube,  f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
 caldb,    s, a, "prod2",,, "Calibration database"
 irf,      s, a, "South_50h",,, "Instrument response function"
+edisp,    b, h, no,,, "Apply energy dispersion?"
 #
 # Upper limit calculation parameters