From b0de6b6f6f083d84b368f9144cd9af3f41778bcc Mon Sep 17 00:00:00 2001 From: mimayer Date: Mon, 6 Jul 2015 16:32:20 +0200 Subject: [PATCH] added support of edisp parameter to ctools and cscripts --- cscripts/cshessobs.py | 35 ++++++++++++++++++++++++++++++++++- cscripts/csresmap.py | 5 ++++- cscripts/cssens.py | 8 +++++--- cscripts/csspec.py | 8 ++++++++ src/ctbutterfly/ctbutterfly.cpp | 25 +++++++++++++++++++++++++ src/ctbutterfly/ctbutterfly.hpp | 1 + src/ctbutterfly/ctbutterfly.par | 1 + src/cterror/cterror.cpp | 25 +++++++++++++++++++++++++ src/cterror/cterror.hpp | 1 + src/cterror/cterror.par | 2 ++ src/cttsmap/cttsmap.cpp | 25 +++++++++++++++++++++++++ src/cttsmap/cttsmap.hpp | 1 + src/cttsmap/cttsmap.par | 1 + src/ctulimit/ctulimit.cpp | 26 +++++++++++++++++++++++++- src/ctulimit/ctulimit.hpp | 1 + src/ctulimit/ctulimit.par | 1 + 16 files changed, 160 insertions(+), 6 deletions(-) diff --git a/cscripts/cshessobs.py b/cscripts/cshessobs.py index 09bf42c..13bf7e0 100755 --- a/cscripts/cshessobs.py +++ b/cscripts/cshessobs.py @@ -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: diff --git a/cscripts/csresmap.py b/cscripts/csresmap.py index 0c061ec..f088a84 100755 --- a/cscripts/csresmap.py +++ b/cscripts/csresmap.py @@ -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 diff --git a/cscripts/cssens.py b/cscripts/cssens.py index 7d8d357..c916953 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -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() diff --git a/cscripts/csspec.py b/cscripts/csspec.py index 41dd06a..2b1a92b 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -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 diff --git a/src/ctbutterfly/ctbutterfly.cpp b/src/ctbutterfly/ctbutterfly.cpp index 396aaa3..21aaf21 100644 --- a/src/ctbutterfly/ctbutterfly.cpp +++ b/src/ctbutterfly/ctbutterfly.cpp @@ -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 save_edisp; + save_edisp.assign(m_obs.size(), false); + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* obs = dynamic_cast(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(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& 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()) { diff --git a/src/ctbutterfly/ctbutterfly.hpp b/src/ctbutterfly/ctbutterfly.hpp index 4011a94..c27ee6d 100644 --- a/src/ctbutterfly/ctbutterfly.hpp +++ b/src/ctbutterfly/ctbutterfly.hpp @@ -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 diff --git a/src/ctbutterfly/ctbutterfly.par b/src/ctbutterfly/ctbutterfly.par index 7432812..6fd17c1 100644 --- a/src/ctbutterfly/ctbutterfly.par +++ b/src/ctbutterfly/ctbutterfly.par @@ -27,6 +27,7 @@ caldb, s, a, "prod2",,, "Calibration database" 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 diff --git a/src/cterror/cterror.cpp b/src/cterror/cterror.cpp index f06f82e..15f0e6d 100644 --- a/src/cterror/cterror.cpp +++ b/src/cterror/cterror.cpp @@ -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 save_edisp; + save_edisp.assign(m_obs.size(), false); + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* obs = dynamic_cast(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(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& 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; diff --git a/src/cterror/cterror.hpp b/src/cterror/cterror.hpp index 206231b..eac6e14 100644 --- a/src/cterror/cterror.hpp +++ b/src/cterror/cterror.hpp @@ -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 diff --git a/src/cterror/cterror.par b/src/cterror/cterror.par index 7f23c68..ea6e626 100644 --- a/src/cterror/cterror.par +++ b/src/cterror/cterror.par @@ -1,3 +1,4 @@ + ############################################################################# # # # CTA Science Analysis Tools # @@ -26,6 +27,7 @@ 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, "prod2",,, "Calibration database" irf, s, a, "South_50h",,, "Instrument response function" +edisp, b, h, no,,, "Apply energy dispersion?" # # Error calculation parameters diff --git a/src/cttsmap/cttsmap.cpp b/src/cttsmap/cttsmap.cpp index a653973..460ec20 100644 --- a/src/cttsmap/cttsmap.cpp +++ b/src/cttsmap/cttsmap.cpp @@ -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 save_edisp; + save_edisp.assign(m_obs.size(), false); + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* obs = dynamic_cast(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(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& 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(); diff --git a/src/cttsmap/cttsmap.hpp b/src/cttsmap/cttsmap.hpp index 89da838..9cf2a8f 100644 --- a/src/cttsmap/cttsmap.hpp +++ b/src/cttsmap/cttsmap.hpp @@ -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 diff --git a/src/cttsmap/cttsmap.par b/src/cttsmap/cttsmap.par index 9343ea4..f3f35f9 100644 --- a/src/cttsmap/cttsmap.par +++ b/src/cttsmap/cttsmap.par @@ -26,6 +26,7 @@ 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" outmap, f, a, "tsmap.fits",,, "Output Test Statistic map" +edisp, b, h, no,,, "Apply energy dispersion?" # # Spatial binning parameters diff --git a/src/ctulimit/ctulimit.cpp b/src/ctulimit/ctulimit.cpp index eb0258a..a858edc 100644 --- a/src/ctulimit/ctulimit.cpp +++ b/src/ctulimit/ctulimit.cpp @@ -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 save_edisp; + save_edisp.assign(m_obs.size(), false); + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* obs = dynamic_cast(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(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& 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; diff --git a/src/ctulimit/ctulimit.hpp b/src/ctulimit/ctulimit.hpp index 2635d2d..155dea1 100644 --- a/src/ctulimit/ctulimit.hpp +++ b/src/ctulimit/ctulimit.hpp @@ -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 diff --git a/src/ctulimit/ctulimit.par b/src/ctulimit/ctulimit.par index e135d2c..88675cf 100644 --- a/src/ctulimit/ctulimit.par +++ b/src/ctulimit/ctulimit.par @@ -25,6 +25,7 @@ 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, "prod2",,, "Calibration database" irf, s, a, "South_50h",,, "Instrument response function" +edisp, b, h, no,,, "Apply energy dispersion?" # # Upper limit calculation parameters -- libgit2 0.21.4