Merge branch '1290-SpecpointsAndTS' into integration

Jürgen Knödlseder
2 parents 96ec743c 337b228b
Showing 2 changed files with 215 additions and 1 deletions
examples/make_spectrum.py
pyext/ctools/obsutils.py
+#! /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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+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
@@ -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