From 0d13fc6058437b50c52bc7bd35934e43aef26135 Mon Sep 17 00:00:00 2001 From: Rolf Buehler Date: Thu, 10 Jul 2014 15:00:29 +0200 Subject: [PATCH] Spectrum points --- scripts/obsutils.py | 91 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- src/ctlike/ctlike.cpp | 33 +++++++++++++++++++++++++++++++++ src/ctselect/ctselect.cpp | 62 ++++++++++++++++++++++++++++++++++---------------------------- 3 files changed, 157 insertions(+), 29 deletions(-) diff --git a/scripts/obsutils.py b/scripts/obsutils.py index 96be22d..00809fc 100755 --- a/scripts/obsutils.py +++ b/scripts/obsutils.py @@ -20,7 +20,7 @@ # ========================================================================== import ctools import gammalib - +import numpy # ===================== # # Simulate observations # @@ -516,3 +516,92 @@ 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 aCTA 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 +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 diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp index f9e3dff..1060f6a 100644 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -267,6 +267,39 @@ void ctlike::run(void) // Optimize model parameters using LM optimizer optimize_lm(); + + // Store original models and max. likelihood + double logL_src=m_logL; + GModels* models_orig = m_obs.models().clone(); + + //Store models whith free parameters + std::vector free_srcs; + GModels models = m_obs.models(); + for (int ii=0; ii < models.size();ii++){ + GModel* model = models[ii]; + for(int jj=0; jj < model->size();jj++){ + GModelPar par = model->at(jj); + if (par.is_free()){ + free_srcs.push_back(model->name()); + break; + } + } + } + //Loop over stored models, remove source and refit + for (int ii=0; ii < free_srcs.size();ii++){ + models.remove(free_srcs[ii]); + m_obs.models(models); + optimize_lm(); + double logL_nosrc=m_logL; + double TS = 2* (logL_src-logL_nosrc); + std::string name = free_srcs[ii]; + (*models_orig)[name]->ts(TS); + models = *models_orig; + } + // Restore best fit values + // TODO: store fit results before and avoid refitting + m_obs.models(*models_orig); + optimize_lm(); // Compute number of observed events in all observations double num_events = 0.0; diff --git a/src/ctselect/ctselect.cpp b/src/ctselect/ctselect.cpp index b01b300..341fd81 100644 --- a/src/ctselect/ctselect.cpp +++ b/src/ctselect/ctselect.cpp @@ -511,7 +511,8 @@ void ctselect::select_events(GCTAObservation* obs, const std::string& filename) // Set requested selections bool select_time = (m_tmin != 0.0 || m_tmax != 0.0); - + bool select_roi = ((m_dec != -1.0 && m_ra != -1.0 && m_rad != -1.0) || m_usepnt); + // Set RA/DEC selection double ra = m_ra; double dec = m_dec; @@ -520,7 +521,7 @@ void ctselect::select_events(GCTAObservation* obs, const std::string& filename) ra = pnt.dir().ra_deg(); dec = pnt.dir().dec_deg(); } - + // Set time selection interval. We make sure here that the time selection // interval cannot be wider than the GTIs covering the data. This is done // using GGti's reduce() method. @@ -565,26 +566,29 @@ void ctselect::select_events(GCTAObservation* obs, const std::string& filename) log << gammalib::parformat("Energy range"); log << m_emin << " - " << m_emax << " TeV" << std::endl; } - if (selection.length() > 0) { - selection += " && "; - } // Make ROI selection - sprintf(cra, "%.6f", ra); - sprintf(cdec, "%.6f", dec); - sprintf(crad, "%.6f", m_rad); - selection += "ANGSEP("+std::string(cra)+"," + - std::string(cdec)+",RA,DEC) <= " + - std::string(crad); - if (logTerse()) { - log << gammalib::parformat("Acceptance cone centre"); - log << "RA=" << ra << ", DEC=" << dec << " deg" << std::endl; - log << gammalib::parformat("Acceptance cone radius"); - log << m_rad << " deg" << std::endl; - } - if (logTerse()) { - log << gammalib::parformat("cfitsio selection"); - log << selection << std::endl; + if (select_roi){ + if (selection.length() > 0) { + selection += " && "; + } + + sprintf(cra, "%.6f", ra); + sprintf(cdec, "%.6f", dec); + sprintf(crad, "%.6f", m_rad); + selection += "ANGSEP("+std::string(cra)+"," + + std::string(cdec)+",RA,DEC) <= " + + std::string(crad); + if (logTerse()) { + log << gammalib::parformat("Acceptance cone centre"); + log << "RA=" << ra << ", DEC=" << dec << " deg" << std::endl; + log << gammalib::parformat("Acceptance cone radius"); + log << m_rad << " deg" << std::endl; + } + if (logTerse()) { + log << gammalib::parformat("cfitsio selection"); + log << selection << std::endl; + } } // Add additional expression @@ -645,11 +649,13 @@ void ctselect::select_events(GCTAObservation* obs, const std::string& filename) // Get CTA event list pointer GCTAEventList* list = static_cast(const_cast(obs->events())); - - // Set ROI - GCTAInstDir instdir; - instdir.dir().radec_deg(ra, dec); - list->roi(GCTARoi(instdir, m_rad)); + + if (select_roi){ + // Set ROI + GCTAInstDir instdir; + instdir.dir().radec_deg(ra, dec); + list->roi(GCTARoi(instdir, m_rad)); + } // Set GTI list->gti(gti); @@ -689,9 +695,9 @@ void ctselect::init_members(void) m_outfile.clear(); m_prefix.clear(); m_usepnt = false; - m_ra = 0.0; - m_dec = 0.0; - m_rad = 0.0; + m_ra = -1.0; + m_dec = -1.0; + m_rad = -1.0; m_tmin = 0.0; m_tmax = 0.0; m_emin = 0.0; -- libgit2 0.21.4