diff --git a/pyext/ctools/obsutils.py b/pyext/ctools/obsutils.py index a15f01d..96be22d 100755 --- a/pyext/ctools/obsutils.py +++ b/pyext/ctools/obsutils.py @@ -516,97 +516,3 @@ 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=numpy.log10(emin) - lgemax=numpy.log10(emax) - lgbins=[] - lgedges=numpy.arange(lgemin,lgemax,(lgemax-lgemin)/nbins) - lgedges=numpy.append(lgedges,[lgemax]) - for ii in range(len(lgedges)-1): - lgbins.append([lgedges[ii],lgedges[ii+1]]) - ebins=numpy.power(10,lgbins) - return ebins - -#Please note that specpoints can make spectra for any given -#spectral model, as long as it has a "Prefactor" parameter -#which scales the funtion -def specpoints(obs,source,ebins): - """ - Returns spectral points in the given binning for - the passed observations and source. - - 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],..] - Keywords: - """ - #Fix all parameters except prefactor - for model in obs.models(): - for par in range(model.size()): - model[par].fix() - - opt_spectral = obs.models()[source].spectral() - opt_spectral["Prefactor"].free() - - #Define spectrum dictionary, units are in format readable by astropy - pars = {"value":[],"eu_value":[],"ed_value":[],"unit":[]} - spec = {"energy":pars.copy(),"flux":pars.copy(),"TS":[]} - - #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 = numpy.power(10,(numpy.log10(ebin[1])+numpy.log10(ebin[0]))/2.) - allflux = opt_spectral.eval(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(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.run() - obsclone = select.obs() - - #Refit and write out to spectrum - binlike=fit(obsclone,log=False, debug=False, edisp=False) - binopt_spectral = binlike.obs().models()[source].spectral() - binflux=binopt_spectral.eval(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(0)) - if binopt_spectral["Prefactor"].value()>0: - relerr=binopt_spectral["Prefactor"].error()/binopt_spectral["Prefactor"].value() - else: - relerr=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"].append("1e12 eV") - - 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"].append("cm-2 1e-12 eV-1") - - spec["TS"].append(binlike.obs().models()[source].ts()) - - #Return spectrum - return spec