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 32504b4..7b6bcba 100755 --- a/pyext/ctools/obsutils.py +++ b/pyext/ctools/obsutils.py @@ -522,109 +522,94 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): # Return observation definition list return obsdeflist - -# ================================================= # -# Make spectrum for one Source in a CTA Observation # -# ================================================= # -def getlogbins(nbins=12,emin=0.1,emax=100): - """ - Returns list with logarithmiscenergy bins in - format [[emin1,emax1],[emin2,emax2],..] - - Parameters: - nbins Number of bins - emin Minimum energy - emax Maximum energy - """ - lgemin=math.log10(emin) - lgemax=math.log10(emax) - lgsize=(lgemax-lgemin)/nbins - ebins=[] - for ii in range(0,nbins): - ebins.append([math.pow(10,lgemin),math.pow(10,lgemin+lgsize)]) - lgemin=lgemin+lgsize - return ebins - -# 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 -def specpoints(obs,source,ebins): + +# ================================================================= # +# 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. + 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 - Pointing direction - source - Source name for which spectrum should be done - ebins - Binning passed as a list in units of TeV as - [[emin1,emax1],[emin2,emax2],..] + 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" + # 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'" + print("Error: your spectral model of "+source+"does not"+\ + "have a scaling parameter named 'Prefactor'"+\ + " or 'Integral'") return - #Fix all parameters except prefactor + # 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 + # 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 ebin in ebins: - #Get logarithmic bin center and flux value of all energy fit - ecenter = math.pow(10,(math.log10(ebin[1])+math.log10(ebin[0]))/2.) - allflux = opt_spectral.eval(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(0)) + # Loop over energy bins, fit and write result to spectrum + for i in range(len(ebounds)): - #Clone observations and reset energy thresholds + # 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) - select["tmax"].real(0) - select["emin"].real(ebin[0]) - select["emax"].real(ebin[1]) + 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() - obsclone = select.obs() - #Refit and write out to spectrum - binlike=fit(obsclone,log=False, debug=False, edisp=False,tscalc=True) + # 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(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(0)) - if binopt_spectral[normpar].value()>0: - relerr=binopt_spectral[normpar].error()/binopt_spectral[normpar].value() + 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 + relerr = 0.0 - #Store and convert to TeV from MeV - spec["energy"]["value"].append(ecenter) - spec["energy"]["ed_value"].append(ecenter-ebin[0]) - spec["energy"]["eu_value"].append(ebin[1]-ecenter) - spec["energy"]["unit"]="1e12 eV" + # 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*1e6) - spec["flux"]["ed_value"].append(relerr*binflux*1e6) - spec["flux"]["eu_value"].append(relerr*binflux*1e6) - spec["flux"]["unit"]="cm-2 1e-12 eV-1" + 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()) - print spec["energy"]["value"] - print spec["flux"]["value"] - #Return spectrum + # Debugging??? + #print(spec["energy"]["value"]) + #print(spec["flux"]["value"]) + + # Return spectrum return spec