Spectrum points

Rolf Buehler
1 parent 0d9c23cf
Showing 3 changed files with 157 additions and 29 deletions
scripts/obsutils.py
src/ctlike/ctlike.cpp
src/ctselect/ctselect.cpp
@@ -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
@@ -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<std::string> 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;
@@ -511,7 +511,8 @@ void ctselect::select_events(GCTAObservation* obs, const std::string&amp; 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&amp; 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&amp; 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&amp; filename)
     // Get CTA event list pointer
     GCTAEventList* list =
         static_cast<GCTAEventList*>(const_cast<GEvents*>(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;