From e3a0f9eb0ebfeaf22e0e63006ae0a32189d6ece7 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 2 Mar 2015 21:11:33 +0100 Subject: [PATCH] Implement computation of integral sensitivity --- doc/source/reference_manual/cssens.rst | 9 +++------ examples/show_sensitivity.py | 49 +++++++++++++++++++++++++++++++++++++------------ scripts/cssens.py | 69 +++++++++++++++++++++++++++++++++++++++++++++------------------------ 3 files changed, 85 insertions(+), 42 deletions(-) mode change 100644 => 100755 examples/show_sensitivity.py diff --git a/doc/source/reference_manual/cssens.rst b/doc/source/reference_manual/cssens.rst index dd53ae2..36f4033 100644 --- a/doc/source/reference_manual/cssens.rst +++ b/doc/source/reference_manual/cssens.rst @@ -45,12 +45,6 @@ General parameters ``outfile [file]`` ASCII file containing the sensitivity values. -``type [string]`` - Source model type. - -``offset [real]`` - Source offset angle (in degrees). - ``duration [real]`` Effective exposure time (in seconds). @@ -75,6 +69,9 @@ General parameters ``(binsz = 0.05) [real]`` Pixel size for binned analysis. +``(type = Differential) [string]`` + Sensitivity type. + ``(sigma = 5.0) [real]`` Significance threshold. diff --git a/examples/show_sensitivity.py b/examples/show_sensitivity.py old mode 100644 new mode 100755 index 1cd538a..4ca60bb --- a/examples/show_sensitivity.py +++ b/examples/show_sensitivity.py @@ -46,13 +46,35 @@ def read_sensitivity(filename): values.append(float(csv[row+1,column])) sensitivity[name] = values + # Check where file contains differential or integral sensitivity + mode = "Integral" + if sensitivity.has_key("emax"): + emax_ref = -1.0 + for value in sensitivity["emax"]: + if emax_ref < 0.0: + emax_ref = value + elif emax_ref != value: + mode = "Differential" + break + + # Add mode + sensitivity["mode"] = mode + # Add linear energy values - if sensitivity.has_key("loge"): - name = 'energy' - values = [] - for value in sensitivity["loge"]: - values.append(math.pow(10.0, value)) - sensitivity[name] = values + if mode == "Differential": + if sensitivity.has_key("loge"): + name = 'energy' + values = [] + for value in sensitivity["loge"]: + values.append(math.pow(10.0, value)) + sensitivity[name] = values + else: + if sensitivity.has_key("emin"): + name = 'energy' + values = [] + for value in sensitivity["emin"]: + values.append(value) + sensitivity[name] = values # Return return sensitivity @@ -61,10 +83,13 @@ def read_sensitivity(filename): # ===================== # # Show sensitivity data # # ===================== # -def show_sensitivity(sensitivity, title="Sensitivity"): +def show_sensitivity(sensitivity, filename): """ Plot sensitivity. """ + # Build title + title = sensitivity['mode'] + ' sensitivity ('+filename+')' + # Create figure plt.figure() @@ -78,7 +103,10 @@ def show_sensitivity(sensitivity, title="Sensitivity"): # Set labels plt.xlabel("Energy (TeV)") - plt.ylabel(r"Sensitivity (erg cm$^{-2}$ s$^{-1}$)") + if sensitivity['mode'] == "Differential": + plt.ylabel(r"E$\times$ F(E) (erg cm$^{-2}$ s$^{-1}$)") + else: + plt.ylabel(r"E$\times$ F($>$E) (erg cm$^{-2}$ s$^{-1}$)") # Show plot plt.show() @@ -105,11 +133,8 @@ if __name__ == '__main__': else: filename = sys.argv[1] - # Set title - title = 'Differential sensitivity ('+filename+')' - # Read sensitivity data sensitivity = read_sensitivity(filename) # Show sensitivity data - show_sensitivity(sensitivity, title=title) + show_sensitivity(sensitivity, filename) diff --git a/scripts/cssens.py b/scripts/cssens.py index 8fb2420..2103148 100755 --- a/scripts/cssens.py +++ b/scripts/cssens.py @@ -47,9 +47,9 @@ class cssens(ctools.cscript): self.version = "1.0.0" # Initialise some parameters - self.obs = None - self.m_ra = None - self.m_dec = None + self.obs = None + self.m_ra = None + self.m_dec = None # Make sure that parfile exists file = self.parfile() @@ -113,6 +113,7 @@ class cssens(ctools.cscript): pars.append(gammalib.GApplicationPar("enumbins","i","h","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("type","s","h","Differential","Differential|Integral","","Sensitivity type")) pars.append(gammalib.GApplicationPar("sigma","r","h","5.0","","","Significance threshold")) pars.append(gammalib.GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations")) pars.append(gammalib.GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average")) @@ -147,17 +148,18 @@ class cssens(ctools.cscript): # Read parameters for binned if requested self.m_enumbins = self["enumbins"].integer() if not self.m_enumbins == 0: - self.m_npix = self["npix"].integer() - self.m_binsz = self["binsz"].real() + self.m_npix = self["npix"].integer() + self.m_binsz = self["binsz"].real() else: # Set dummy values (required by obsutils) - self.m_npix = 0 + self.m_npix = 0 self.m_binsz = 0.0 # Read remaining parameters 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() + self.m_type = self["type"].string() # Set some fixed parameters self.m_log = False # Logging in client tools @@ -219,20 +221,28 @@ class cssens(ctools.cscript): if self.logTerse(): self.log("\n") self.log.header1("Sensitivity determination") + self.log.parformat("Type") + self.log(self.m_type) + self.log("\n") # Loop over energy bins for ieng in range(self.m_ebounds.size()): # Set energies - emean = self.m_ebounds.elogmean(ieng).TeV() - emin = self.m_ebounds.emin(ieng).TeV() - emax = self.m_ebounds.emax(ieng).TeV() - - # Set energy boundaries - self.set_obs_ebounds(emin, emax) + if self.m_type == "Differential": + emin = self.m_ebounds.emin(ieng) + emax = self.m_ebounds.emax(ieng) + elif self.m_type == "Integral": + emin = self.m_ebounds.emin(ieng) + emax = self.m_ebounds.emax() + else: + msg = "Invalid sensitivity type \""+self.m_type+"\" encountered."+ \ + " Either use \"Differential\" or \"Integral\"." + raise gammalib.GException.invalid_value("cssens", msg) # Determine sensitivity - result = self.get_sensitivity(self.obs, bkg_model, full_model) + result = self.get_sensitivity(self.obs, emin, emax, \ + bkg_model, full_model) # Write results if ieng == 0: @@ -305,13 +315,16 @@ class cssens(ctools.cscript): def set_obs_ebounds(self, emin, emax): """ Set energy boundaries for observation in container. + + Parameters: + emin - Minimum energy + emax - Maximum energy """ # Loop over all observations in container for obs in self.obs: # Set energy boundaries - ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \ - gammalib.GEnergy(emax, "TeV")) + ebounds = gammalib.GEbounds(emin, emax) obs.events().ebounds(ebounds) # Return @@ -382,27 +395,35 @@ class cssens(ctools.cscript): # Return models return full_model, bkg_model - def get_sensitivity(self, obs, bkg_model, full_model): + def get_sensitivity(self, obs, emin, emax, bkg_model, full_model): """ - Determine sensitivity for a given observations. + Determine sensitivity for given observations. Parameters: obs - Observation container + emin - Minimum energy for fitting and flux computation + emax - Maximum energy for fitting and flux computation bkg_model - Background model full_model - Source model """ # Set TeV->erg conversion factor tev2erg = 1.6021764 + + # Set energy boundaries + self.set_obs_ebounds(emin, emax) # Determine energy boundaries from first observation in # the container - for run in obs: - emin = run.events().ebounds().emin() - emax = run.events().ebounds().emax() - loge = math.log10(math.sqrt(emin.TeV()*emax.TeV())) - erg_mean = math.pow(10.0, loge) * tev2erg - erg_width = (emax.TeV()-emin.TeV()) * tev2erg - break + #for run in obs: + # emin = run.events().ebounds().emin() + # emax = run.events().ebounds().emax() + # loge = math.log10(math.sqrt(emin.TeV()*emax.TeV())) + # erg_mean = math.pow(10.0, loge) * tev2erg + # erg_width = (emax.TeV()-emin.TeV()) * tev2erg + # break + loge = math.log10(math.sqrt(emin.TeV()*emax.TeV())) + erg_mean = math.pow(10.0, loge) * tev2erg + erg_width = (emax.TeV()-emin.TeV()) * tev2erg # Compute Crab unit (this is the factor with which the Prefactor needs # to be multiplied to get 1 Crab -- libgit2 0.21.4