diff --git a/pyext/ctools/obsutils.py b/pyext/ctools/obsutils.py
index 6f632c5..ca125c8 100755
--- a/pyext/ctools/obsutils.py
+++ b/pyext/ctools/obsutils.py
@@ -571,90 +571,3 @@ 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"] = "UNDEF"
-        select["dec"] = "UNDEF"
-        select["rad"] = "UNDEF"
-        select["tmin"] = 0
-        select["tmax"] = 0
-        select["emin"] = ebounds.emin(i).TeV()
-        select["emax"] = ebounds.emax(i).TeV()
-        select.run()
-        
-        # Refit and write out to spectrum
-        binlike         = fit(select.obs(), log=False, debug=False, edisp=False)
-        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())
-    
-    # Return spectrum
-    return spec