diff --git a/examples/show_sensitivity.py b/examples/show_sensitivity.py index 4ca60bb..28c6ccf 100755 --- a/examples/show_sensitivity.py +++ b/examples/show_sensitivity.py @@ -99,7 +99,7 @@ def show_sensitivity(sensitivity, filename): plt.grid() # Show differential sensitivity - plt.plot(sensitivity['energy'], sensitivity['eflux'], 'ro-') + plt.plot(sensitivity['energy'], sensitivity['sensitivity'], 'ro-') # Set labels plt.xlabel("Energy (TeV)") diff --git a/pyext/ctools/obsutils.py b/pyext/ctools/obsutils.py index 2ad0d87..3224907 100755 --- a/pyext/ctools/obsutils.py +++ b/pyext/ctools/obsutils.py @@ -26,7 +26,7 @@ import math # ===================== # # Simulate observations # # ===================== # -def sim(obs, log=False, debug=False, edisp=False, seed=0, nbins=0, +def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, binsz=0.05, npix=200, proj="TAN", coord="GAL"): """ Simulate events for all observations in the container. @@ -55,6 +55,9 @@ def sim(obs, log=False, debug=False, edisp=False, seed=0, nbins=0, # Optionally switch-on debugging model if debug: sim["debug"].boolean(True) + + # Set chatter level + sim["chatter"].chatter = chatter # Run ctobssim application. This will loop over all observations in the # container and simulation the events for each observation. Note that @@ -101,6 +104,9 @@ def sim(obs, log=False, debug=False, edisp=False, seed=0, nbins=0, if debug: bin["debug"].boolean(True) + # Set chatter level + bin["chatter"].chatter = chatter + # Run ctbin application. This will loop over all observations in # the container and bin the events in counts maps bin.run() @@ -127,7 +133,7 @@ def sim(obs, log=False, debug=False, edisp=False, seed=0, nbins=0, # ================ # # Fit observations # # ================ # -def fit(obs, log=False, debug=False, edisp=False): +def fit(obs, log=False, debug=False, chatter=2, edisp=False): """ Perform maximum likelihood fitting of observations in the container. @@ -148,6 +154,9 @@ def fit(obs, log=False, debug=False, edisp=False): # Optionally switch-on debugging model if debug: like["debug"].boolean(True) + + # Set chatter level + like["chatter"].chatter = chatter # Optionally apply energy dispersion like["edisp"].boolean(edisp) diff --git a/scripts/cssens.py b/scripts/cssens.py index 2103148..f4e8503 100755 --- a/scripts/cssens.py +++ b/scripts/cssens.py @@ -200,7 +200,8 @@ class cssens(ctools.cscript): self.log("\n") # Initialise script - colnames = ['loge', 'emin', 'emax', 'crab', 'pflux', 'eflux', 'diffSens'] + colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux', \ + 'energy_flux', 'sensitivity'] results = [] # Initialise models @@ -412,18 +413,10 @@ class cssens(ctools.cscript): # 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 + # Determine energy boundaries from first observation in the container loge = math.log10(math.sqrt(emin.TeV()*emax.TeV())) - erg_mean = math.pow(10.0, loge) * tev2erg - erg_width = (emax.TeV()-emin.TeV()) * tev2erg + e_mean = math.pow(10.0, loge) + erg_mean = e_mean * tev2erg # Compute Crab unit (this is the factor with which the Prefactor needs # to be multiplied to get 1 Crab @@ -446,12 +439,12 @@ class cssens(ctools.cscript): self.log("\n") # Initialise loop - flux_value = [] - pflux_value = [] - eflux_value = [] - diffSens_value = [] - iter = 0 - test_cflux = 0.1 # This is the initial test flux in Crab units (100 mCrab) + crab_flux_value = [] + photon_flux_value = [] + energy_flux_value = [] + sensitivity_value = [] + iter = 0 + test_crab_flux = 0.1 # Initial test flux in Crab units (100 mCrab) # Loop until we break while True: @@ -463,13 +456,13 @@ class cssens(ctools.cscript): if self.logExplicit(): self.log.header2("Iteration "+str(iter)) - # Set test flux + # Set source model. crab_prefactor is the Prefactor that + # corresponds to 1 Crab src_model = full_model.copy() crab_prefactor = src_model[self.m_srcname]['Prefactor'].value() * crab_unit - src_model[self.m_srcname]['Prefactor'].value(crab_prefactor * test_cflux) + src_model[self.m_srcname]['Prefactor'].value(crab_prefactor * test_crab_flux) obs.models(src_model) - # Simulate events sim = obsutils.sim(obs, nbins=self.m_enumbins, seed=iter, \ binsz=self.m_binsz, npix=self.m_npix, \ @@ -572,48 +565,45 @@ class cssens(ctools.cscript): continue # Get fitted Crab, photon and energy fluxes - cflux = result_all[self.m_srcname]['Prefactor'].value() / crab_prefactor - cflux_err = result_all[self.m_srcname]['Prefactor'].error() / crab_prefactor - pflux = result_all[self.m_srcname].spectral().flux(emin, emax) - eflux = result_all[self.m_srcname].spectral().eflux(emin, emax) + crab_flux = result_all[self.m_srcname]['Prefactor'].value() / crab_prefactor + crab_flux_err = result_all[self.m_srcname]['Prefactor'].error() / crab_prefactor + photon_flux = result_all[self.m_srcname].spectral().flux(emin, emax) + energy_flux = result_all[self.m_srcname].spectral().eflux(emin, emax) - # Compute differential sensitivity - diffSens = pflux / erg_width * erg_mean*erg_mean + # Compute differential sensitivity in unit erg/cm2/s + energy = gammalib.GEnergy(e_mean, "TeV") + time = gammalib.GTime() + sensitivity = result_all[self.m_srcname].spectral().eval(energy, time) * \ + erg_mean*erg_mean * 1.0e6 # Compute flux correction factor based on average TS correct = 1.0 if ts > 0: correct = math.sqrt(self.m_ts_thres/ts) - #if self.m_use_ts: - # if ts > 0: - # correct = math.sqrt(self.m_ts_thres/ts) - #else: - # if cflux > 0 and cflux_err > 0: - # correct = math.sqrt(self.m_ts_thres)/(cflux/cflux_err) # Compute extrapolated fluxes - flux = correct * cflux - pflux = correct * pflux - eflux = correct * eflux - diffSens = correct * diffSens - flux_value.append(flux) - pflux_value.append(pflux) - eflux_value.append(eflux) - diffSens_value.append(diffSens) + crab_flux = correct * crab_flux + photon_flux = correct * photon_flux + energy_flux = correct * energy_flux + sensitivity = correct * sensitivity + crab_flux_value.append(crab_flux) + photon_flux_value.append(photon_flux) + energy_flux_value.append(energy_flux) + sensitivity_value.append(sensitivity) # Write background and test source fit results if self.logExplicit(): self.log.parformat("Photon flux") - self.log(pflux) + self.log(photon_flux) self.log(" ph/cm2/s\n") self.log.parformat("Energy flux") - self.log(eflux) + self.log(energy_flux) self.log(" erg/cm2/s\n") self.log.parformat("Crab flux") - self.log(cflux*1000.0) + self.log(crab_flux*1000.0) self.log(" mCrab\n") self.log.parformat("Differential sensitivity") - self.log(diffSens) + self.log(sensitivity) self.log(" erg/cm2/s\n") for model in result_all: self.log.parformat("Model") @@ -629,40 +619,40 @@ class cssens(ctools.cscript): self.log("corr=") self.log(correct) self.log(" ") - self.log(pflux) + self.log(photon_flux) self.log(" ph/cm2/s = ") - self.log(eflux) + self.log(energy_flux) self.log(" erg/cm2/s = ") - self.log(cflux*1000.0) + self.log(crab_flux*1000.0) self.log(" mCrab = ") - self.log(diffSens) + self.log(sensitivity) self.log(" erg/cm2/s\n") # Compute sliding average of extrapolated fitted prefactor, # photon and energy flux. This damps out fluctuations and # improves convergence - flux = 0.0 - pflux = 0.0 - eflux = 0.0 - diffSens = 0.0 - num = 0.0 + crab_flux = 0.0 + photon_flux = 0.0 + energy_flux = 0.0 + sensitivity = 0.0 + num = 0.0 for k in range(self.m_num_avg): - inx = len(flux_value) - k - 1 + inx = len(crab_flux_value) - k - 1 if inx >= 0: - flux += flux_value[inx] - pflux += pflux_value[inx] - eflux += eflux_value[inx] - diffSens += diffSens_value[inx] + crab_flux += crab_flux_value[inx] + photon_flux += photon_flux_value[inx] + energy_flux += energy_flux_value[inx] + sensitivity += sensitivity_value[inx] num += 1.0 - flux /= num - pflux /= num - eflux /= num - diffSens /= num + crab_flux /= num + photon_flux /= num + energy_flux /= num + sensitivity /= num # Compare average flux to last average if iter > self.m_num_avg: - if test_cflux > 0: - ratio = flux/test_cflux + if test_crab_flux > 0: + ratio = crab_flux/test_crab_flux # We have 2 convergence criteria: # 1. The average flux does not change @@ -678,7 +668,7 @@ class cssens(ctools.cscript): break # Use average for next iteration - test_cflux = flux + test_crab_flux = crab_flux # Exit loop if number of trials exhausted if (iter >= self.m_max_iter): @@ -693,16 +683,16 @@ class cssens(ctools.cscript): self.log(ts) self.log("\n") self.log.parformat("Photon flux") - self.log(pflux) + self.log(photon_flux) self.log(" ph/cm2/s\n") self.log.parformat("Energy flux") - self.log(eflux) + self.log(energy_flux) self.log(" erg/cm2/s\n") self.log.parformat("Crab flux") - self.log(cflux*1000.0) + self.log(crab_flux*1000.0) self.log(" mCrab\n") self.log.parformat("Differential sensitivity") - self.log(diffSens) + self.log(sensitivity) self.log(" erg/cm2/s\n") self.log.parformat("Number of simulated events") self.log(nevents) @@ -736,8 +726,9 @@ class cssens(ctools.cscript): # Store result result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \ - 'crab': flux, 'pflux': pflux, 'eflux': eflux, \ - 'diffSens': diffSens} + 'crab_flux': crab_flux, 'photon_flux': photon_flux, \ + 'energy_flux': energy_flux, \ + 'sensitivity': sensitivity} # Return result return result