diff --git a/examples/make_spectrum.py b/examples/make_spectrum.py new file mode 100755 index 0000000..83e8577 --- /dev/null +++ b/examples/make_spectrum.py @@ -0,0 +1,117 @@ +#! /usr/bin/env python +# ========================================================================== +# This script shows how to make a spectrum using obsutils. +# +# Copyright (C) 2014 Juergen Knoedlseder +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import gammalib +import ctools +try: + import matplotlib.pyplot as plt + has_matplotlib = True +except: + has_matplotlib = False + + +# =============== # +# Make a spectrum # +# =============== # +def make_spectrum(): + """ + """ + # Set script parameters + model_name = "${CTOOLS}/share/models/crab.xml" + caldb = "dummy" + irf = "cta_dummy_irf" + ra = 83.63 + dec = 22.01 + rad_sim = 3.0 + tstart = 0.0 + tstop = 1800.0 + emin = 0.1 + emax = 100.0 + + # Simulate events + sim = ctools.ctobssim() + sim["infile"].filename(model_name) + sim["caldb"].string(caldb) + sim["irf"].string(irf) + sim["ra"].real(ra) + sim["dec"].real(dec) + sim["rad"].real(rad_sim) + sim["tmin"].real(tstart) + sim["tmax"].real(tstop) + sim["emin"].real(emin) + sim["emax"].real(emax) + sim.run() + + # Generate an energy binning + e_min = gammalib.GEnergy(emin, "TeV") + e_max = gammalib.GEnergy(emax, "TeV") + ebounds = gammalib.GEbounds(10, e_min, e_max) + + # Generate spectral points + spectrum = ctools.obsutils.spectrum(sim.obs(), "Crab", ebounds) + + # Return spectrum + return spectrum + + +# ============= # +# Plot spectrum # +# ============= # +def plot_spectrum(spectrum): + """ + Plot spectrum. + """ + # Create figure + plt.figure(1) + plt.title("Crab spectrum") + + # Plot spectrum + plt.loglog(spectrum['energy']['value'], \ + spectrum['flux']['value'], 'ro', label='Crab') + plt.errorbar(spectrum['energy']['value'], \ + spectrum['flux']['value'], \ + spectrum['flux']['ed_value'], fmt=None, ecolor='r') + + # Put labels + plt.xlabel("Energy ("+spectrum['energy']['unit']+")") + plt.ylabel("Flux ("+spectrum['flux']['unit']+")") + plt.legend(loc="lower left") + + # Show spectrum + plt.show() + + # Return + return + + +#==========================# +# Main routine entry point # +#==========================# +if __name__ == '__main__': + """ + Generate a Crab spectrum. + """ + # Generate spectrum + spectrum = make_spectrum() + + # Plot spectrum + if has_matplotlib: + plot_spectrum(spectrum) + \ No newline at end of file diff --git a/pyext/ctools/obsutils.py b/pyext/ctools/obsutils.py index 96be22d..7b6bcba 100755 --- a/pyext/ctools/obsutils.py +++ b/pyext/ctools/obsutils.py @@ -20,6 +20,7 @@ # ========================================================================== import ctools import gammalib +import math # ===================== # @@ -124,7 +125,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, edisp=False, tscalc=False): """ Perform maximum likelihood fitting of observations in the container. @@ -145,9 +146,14 @@ def fit(obs, log=False, debug=False, edisp=False): # Optionally switch-on debugging model if debug: like["debug"].boolean(True) + + # Optionally turn on TS calculation + if tscalc: + like["tscalc"].boolean(True) # Optionally apply energy dispersion like["edisp"].boolean(edisp) + like["edisp"].boolean(edisp) # Run ctlike application. like.run() @@ -516,3 +522,94 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): # Return observation definition list return obsdeflist + + +# ================================================================= # +# Generate a spectrum by fitting an observation in energy intervals # +# ================================================================= # +def spectrum(obs, source, ebounds): + """ + Returns spectral points in the given binning for the passed observations + and source. + + Please note that specpoints can make spectra for any given spectral + model, as long as it has a "Prefactor" or "Integral" parameter which + scales the function + + Parameters: + obs - Observation container + source - Source name for which spectrum should be computed + ebounds - GEbounds energy boundaries + Keywords: + None + """ + # Check if parameters Prefactor or Integral exist to scale mdoel + normpar = "Prefactor" + srcmodel = obs.models()[source] + if not srcmodel.has_par(normpar): + if srcmodel.has_par("Integral"): + normpar = "Integral" + else: + print("Error: your spectral model of "+source+"does not"+\ + "have a scaling parameter named 'Prefactor'"+\ + " or 'Integral'") + return + + # Fix all parameters except prefactor + opt_spectral = srcmodel.spectral() + for model in obs.models(): + for par in range(model.size()): + model[par].fix() + opt_spectral[normpar].free() + + # Define spectrum dictionary, units are in format readable by astropy + spec = {"energy":{"value":[],"eu_value":[],"ed_value":[],"unit":"none"}, + "flux" :{"value":[],"eu_value":[],"ed_value":[],"unit":"none"}, + "TS" :{"value":[],"eu_value":[],"ed_value":[],"unit":"none"}} + + # Loop over energy bins, fit and write result to spectrum + for i in range(len(ebounds)): + + # Get logarithmic bin center and flux value of all energy fit + ecenter = ebounds.elogmean(i) + allflux = opt_spectral.eval(ecenter, gammalib.GTime(0.0)) + + # Clone observations and reset energy thresholds + select = ctools.ctselect(obs) + select["ra"].real(-1.0) + select["dec"].real(-1.0) + select["rad"].real(-1.0) + select["tmin"].real(0.0) + select["tmax"].real(0.0) + select["emin"].real(ebounds.emin(i).TeV()) + select["emax"].real(ebounds.emax(i).TeV()) + select.run() + + # Refit and write out to spectrum + binlike = fit(select.obs(), log=False, debug=False, edisp=False, tscalc=True) + binopt_spectral = binlike.obs().models()[source].spectral() + binflux = binopt_spectral.eval(ecenter, gammalib.GTime(0.0)) + if binopt_spectral[normpar].value() > 0.0: + relerr = binopt_spectral[normpar].error()/binopt_spectral[normpar].value() + else: + relerr = 0.0 + + # Store and convert to TeV from MeV + spec["energy"]["value"].append(ecenter.TeV()) + spec["energy"]["ed_value"].append(ecenter.TeV()-ebounds.emin(i).TeV()) + spec["energy"]["eu_value"].append(ebounds.emax(i).TeV()-ecenter.TeV()) + spec["energy"]["unit"] = "TeV" + + spec["flux"]["value"].append(binflux*1.0e6) + spec["flux"]["ed_value"].append(relerr*binflux*1.0e6) + spec["flux"]["eu_value"].append(relerr*binflux*1.0e6) + spec["flux"]["unit"] = "ph/cm2/s/TeV" + + spec["TS"]["value"].append(binlike.obs().models()[source].ts()) + + # Debugging??? + #print(spec["energy"]["value"]) + #print(spec["flux"]["value"]) + + # Return spectrum + return spec