Merge pull request #123 from mireianievas/master

Several (small) changes

enrico/RunGTlike.py

@@ -29,11 +29,10 @@ def GenAnalysisObjects(config, verbose = 1, xmlfile =""):
     mes = Loggin.Message()
     #check is the summed likelihood method should be used and get the
     #Analysis objects (observation and (Un)BinnedAnalysis objects)
-    SummedLike = config['Spectrum']['SummedLike']
     folder = config['out']
 
     # If there are no xml files, create it and print a warning
-    if (glob.glob("%s*.xml" %config['file']['xml'].replace('.xml','')) is []):
+    if len(glob.glob(config['file']['xml'].replace('.xml','*.xml')))==0:
         mes.warning("Xml not found, creating one for the given config %s" %config['file']['xml'])
         XmlMaker(config)
 
@@ -47,7 +46,6 @@ def GenAnalysisObjects(config, verbose = 1, xmlfile =""):
             from enrico.data import fermievtypes
             mes.info("Breaking the analysis in Front/Back events")
             # Set Summed Likelihood to True
-            config['Spectrum']['SummedLike'] = 'yes'
             oldxml = config['file']['xml']
             for k,TYPE in enumerate(["FRONT", "BACK"]):
                 configs[k] = ConfigObj(config)

enrico/config.py

@@ -1,6 +1,10 @@
 """Central place for config file handling"""
 import sys
 from os.path import join
+try:
+    import pyfits
+except:
+    import astropy.io.fits as pyfits
 from enrico.extern.configobj import ConfigObj, flatten_errors
 from enrico.extern.validate import Validator
 from enrico.environ import CONFIG_DIR, DOWNLOAD_DIR, USE_FULLMISSION_SPACECRAFT
@@ -41,6 +45,35 @@ def get_config(infile, configspec=join(CONFIG_DIR, 'default.conf')):
 def get_default_config(configspec=join(CONFIG_DIR, 'default.conf')):
     return ConfigObj(None, configspec=configspec)
 
+def get_times_from_spacecraft(scfile,target=['tmin','tmax']):
+    tset=False
+    tmin = 239557418
+    tmax = 334165418
+    if tset is False:
+        try:
+            sc = pyfits.open(scfile)
+            if 'tmin' in target:
+                tmin = sc[0].header['TSTART']
+            if 'tmax' in target:
+                tmax = sc[0].header['TSTOP']
+        except:
+            raise
+        else:
+            tset=True
+    if tset is False:
+        try:
+            with open(scfile.replace('@','')) as f:
+                if 'tmin' in target:
+                    sc1 = pyfits.open(f.readlines()[0])
+                    tmin = sc1[0].header['TSTART']
+                if 'tmax' in target:
+                    sc2 = pyfits.open(f.readlines()[-1])
+                    tmax = sc2[0].header['TSTOP']
+        except:
+            raise
+        else:
+            tset=True
+    return([tmin,tmax])
 
 def query_config():
     import os
@@ -107,16 +140,17 @@ def query_config():
 
 #    informations about the time
     config['time'] = {}
-    tmin = raw_input('Start time [239557418] : ')
-    if not(tmin=='') :
+    tmin = raw_input('Start time [-1=START] : ')
+    ft2 = config['file']['spacecraft']
+    if not(tmin=='') and float(tmin)>=0 :
       config['time']['tmin'] = tmin
     else :
-      config['time']['tmin'] = '239557418'
-    tmax = raw_input('End time [334165418] : ')
-    if not(tmax=='') :
+      config['time']['tmin'] = get_times_from_spacecraft(ft2,target=['tmin'])[0]
+    tmax = raw_input('End time [-1=END] : ')
+    if not(tmax=='') and float(tmax)>=0 :
       config['time']['tmax'] = tmax
     else :
-      config['time']['tmax'] = '334165418'
+      config['time']['tmax'] = get_times_from_spacecraft(ft2,target=['tmax'])[1]
 
 #    informations about the energy
     config['energy'] = {}

enrico/config/default.conf

@@ -71,8 +71,8 @@ Submit = option('yes', 'no', default='no')
 	phibins = integer(default=0)
 
 [time]
-	tmin = float(default=239557418)
-	tmax = float(default=239557418)
+	tmin = float(default=-1) #239557418)
+	tmax = float(default=-1) #239557418)
 	file = string(default='')
 	type = option('MET', 'MJD', 'JD', default='MET')
 

enrico/gtfunction.py

@@ -11,7 +11,6 @@
 from GtApp import GtApp
 from enrico import utils
 
-
 class Observation:
     # init function of the Observation class.
     # folder : folder where the produced fits files will be stored.

enrico/lightcurve.py

@@ -303,19 +303,19 @@ def _PlotLC(self,folded=False):
 
         # #change the list into np array
         # TS = np.array(TS)
-        Npred = np.array(Npred)
-        Npred_detected = Npred[Npred_detected_indices]
-        Time = np.array(Time)
-        # TimeErr = np.array(TimeErr)
-        Flux = np.array(Flux)
-        FluxErr = np.array(FluxErr)
+        Npred = np.asarray(Npred)
+        Npred_detected = np.asarray(Npred[Npred_detected_indices])
+        Time = np.asarray(Time)
+        TimeErr = np.asarray(TimeErr)
+        Flux = np.asarray(Flux)
+        FluxErr = np.asarray(FluxErr)
         # Index = np.array(Index)
         # IndexErr = np.array(IndexErr)
         # Cutoff = np.array(Cutoff)
         # CutoffErr = np.array(CutoffErr)
-        FluxForNpred = np.array(FluxForNpred)
+        FluxForNpred = np.asarray(FluxForNpred)
         # FluxErrForNpred = np.array(FluxErrForNpred)
-
+        uplim = np.asarray(uplim,dtype=bool)
         #Plots the diagnostic plots is asked
         # Plots are : Npred vs flux
         #             TS vs Time
@@ -353,10 +353,26 @@ def _PlotLC(self,folded=False):
 
             #plot TS vs Time
             plt.figure()
-            plt.ylim(ymin=min(TS)*0.8,ymax=max(TS)*1.2)
             plt.xlabel(r"Time (s)")
             plt.ylabel(r"Test Statistic")
-            plt.errorbar(Time,TS,xerr=TimeErr,fmt='+',color='black',ls='None')
+            plt.errorbar(x=Time,y=TS,xerr=TimeErr,fmt='+',color='black',ls='None')
+            plt.ylim(ymin=min(TS)*0.8,ymax=max(TS)*1.2)
+            plt.xlim(xmin=max(plt.xlim()[0],1.02*min(Time)-0.02*max(Time)),xmax=min(plt.xlim()[1],1.02*max(Time)-0.02*min(Time)))
+
+            # Move the offset to the axis label
+            ax = plt.gca()
+            ax.get_yaxis().get_major_formatter().set_useOffset(False)
+            offset_factor = int(np.mean(np.log10(np.abs(ax.get_ylim()))))
+            if (offset_factor != 0):
+                ax.set_yticklabels([float(round(k,5)) for k in ax.get_yticks()*10**(-offset_factor)])
+                ax.yaxis.set_label_text(ax.yaxis.get_label_text() + r" [${\times 10^{%d}}$]" %offset_factor)
+
+            # Secondary axis with MJD
+            mjdaxis = plt.twiny()
+            mjdaxis.set_xlim([utils.met_to_MJD(k) for k in ax.get_xlim()])
+            mjdaxis.set_xlabel(r"Time (MJD)")
+            plt.tight_layout()
+
             plt.savefig(LcOutPath+"_TS.png", dpi=150, facecolor='w', edgecolor='w',
                     orientation='portrait', papertype=None, format=None,
                     transparent=False, bbox_inches=None, pad_inches=0.1,
@@ -390,12 +406,50 @@ def _PlotLC(self,folded=False):
         # plt.xlim(xmin=xmin,xmax=xmax)
         #plt.errorbar(Time,Flux,xerr=TimeErr,yerr=FluxErr,fmt='o',color='black',ls='None',uplims=uplim)
         plot_errorbar_withuls(Time,TimeErr,TimeErr,Flux,FluxErr,FluxErr,uplim,bblocks=True)
-
+        plt.ylim(ymin=max(plt.ylim()[0],np.percentile(Flux[~uplim],1)*0.1),
+                 ymax=min(plt.ylim()[1],np.percentile(Flux[~uplim],99)*2.0))
+        plt.xlim(xmin=max(plt.xlim()[0],1.02*min(Time)-0.02*max(Time)),
+                 xmax=min(plt.xlim()[1],1.02*max(Time)-0.02*min(Time)))
+
+        # Move the offset to the axis label
+        ax = plt.gca()
+        ax.get_yaxis().get_major_formatter().set_useOffset(False)
+        offset_factor = int(np.mean(np.log10(np.abs(ax.get_ylim()))))
+        if (offset_factor != 0):
+            ax.set_yticklabels([float(round(k,5)) for k in ax.get_yticks()*10**(-offset_factor)])
+            ax.yaxis.set_label_text(ax.yaxis.get_label_text() + r" [${\times 10^{%d}}$]" %offset_factor)
+
+        # Secondary axis with MJD
+        mjdaxis = plt.twiny()
+        mjdaxis.set_xlim([utils.met_to_MJD(k) for k in ax.get_xlim()])
+        mjdaxis.set_xlabel(r"Time (MJD)")
+        plt.tight_layout()
+
         plt.savefig(LcOutPath+"_LC.png", dpi=150, facecolor='w', edgecolor='w',
                 orientation='portrait', papertype=None, format=None,
                 transparent=False, bbox_inches=None, pad_inches=0.1,
                 frameon=None)
 
+        if self.config["LightCurve"]["SpectralIndex"] == 0 :
+            plt.figure()
+            plt.xlabel(r"Time (s)")
+            plt.ylabel(r"${\rm Index}$")
+            Index = np.asarray(Index)
+            IndexErr = np.asarray(IndexErr)
+            uplimIndex = uplim + Index<0.55
+            plot_errorbar_withuls(Time[~uplimIndex],TimeErr[~uplimIndex],TimeErr[~uplimIndex],
+                                  Index[~uplimIndex],IndexErr[~uplimIndex],IndexErr[~uplimIndex],
+                                  uplimIndex[~uplimIndex],bblocks=True)
+
+            plt.ylim(ymin=max(plt.ylim()[0],np.percentile(Index[~uplimIndex],1)*0.1),ymax=min(plt.ylim()[1],np.percentile(Index[~uplimIndex],99)*2.0))
+            plt.xlim(xmin=max(plt.xlim()[0],1.02*min(Time)-0.02*max(Time)),xmax=min(plt.xlim()[1],1.02*max(Time)-0.02*min(Time)))
+
+            plt.savefig(LcOutPath+"_Index.png", dpi=150, facecolor='w', edgecolor='w',
+                orientation='portrait', papertype=None, format=None,
+                transparent=False, bbox_inches=None, pad_inches=0.1,
+                frameon=None)
+
+
         # compute Fvar and probability of being cst
 
         self.info("Flux vs Time: infos")

enrico/plotting.py

@@ -181,10 +181,11 @@ def CountsPlot(self, Parameter):
         plt.xlabel("E (MeV) ")
         plt.ylabel("Counts / bin")
         plt.errorbar(E,obs,xerr=err_E,yerr=obs_err,fmt='o',color="red",ls='None',label="Data")
-        plt.plot(E,other,'--',color="blue",label=Parameter.srcname)
-        plt.plot(E,src,'-',color="green",label="Other Sources")
-        plt.plot(E,total,'-',ls='None',label="All Sources")
+        plt.plot(E,other,ls='dashed',color="blue",label=Parameter.srcname)
+        plt.plot(E,src,ls='solid',color="green",label="Other Sources")
+        plt.plot(E,total,lw=1.5,ls='solid',label="All Sources")
         plt.legend()
+        plt.tight_layout()
         plt.savefig(filebase + "_CountsPlot.png", dpi=150, facecolor='w', edgecolor='w',
             orientation='portrait', papertype=None, format=None,
             transparent=False, bbox_inches=None, pad_inches=0.1,
@@ -212,8 +213,9 @@ def CountsPlot(self, Parameter):
         plt.ylabel("(counts-model)/model")
         plt.errorbar(E,residual,xerr=err_E,yerr=Dres,fmt='o',color="red",ls='None',label="Data")
         zero = np.zeros(2)
-        Ezero = np.array([0, 1e10])
-        plt.plot(Ezero,zero,'-',color='black')
+        Ezero = np.array([1e-5, 1e10])
+        plt.plot(Ezero,zero,lw=1.5,ls='solid',color='black')
+        plt.tight_layout()
         plt.savefig(filebase + "ResPlot.png", dpi=150, facecolor='w', edgecolor='w',
             orientation='portrait', papertype=None, format=None,
             transparent=False, bbox_inches=None, pad_inches=0.1,
@@ -339,11 +341,13 @@ def plot_errorbar_withuls(x,xerrm,xerrp,y,yerrm,yerrp,uplim,bblocks=False):
     yerrp[uplim] = 0
 
     optimal_markersize = (0.5+4./(1.+np.log10(len(y))))
+    optimal_errorlinewidth = (0.2+2./(1.+4.*np.log10(len(y))))
 
     # Plot the significant points
     plt.errorbar(x[~uplim], y[~uplim],
         xerr=[xerrm[~uplim], xerrp[~uplim]],
         yerr=[yerrm[~uplim], yerrp[~uplim]],
+        lw=optimal_errorlinewidth,
         fmt='o',ms=optimal_markersize,capsize=0,zorder=10,
         color='black',ls='None',uplims=False,label='LAT data')
 
@@ -354,19 +358,23 @@ def plot_errorbar_withuls(x,xerrm,xerrp,y,yerrm,yerrp,uplim,bblocks=False):
             xerr=[xerrm[uplim], xerrp[uplim]],
             yerr=[yerrm[uplim], yerrp[uplim]],
             fmt='o',markersize=0,capsize=0,zorder=-1,
+            lw=optimal_errorlinewidth,
             color='0.50',ls='None',lolims=False)
         plt.errorbar(x[uplim], 0.8*y[uplim],
             yerr=[0.2*y[uplim], 0.2*y[uplim]],
             fmt='o',markersize=0,capsize=optimal_markersize/1.5,zorder=-1,
+            lw=optimal_errorlinewidth,
             color='0.50',ls='None',lolims=True)
     else:
         plt.errorbar(x[uplim], y[uplim],
             xerr=[xerrm[uplim], xerrp[uplim]],
             yerr=[yerrm[uplim], yerrp[uplim]],
+            lw=optimal_errorlinewidth,
             fmt='o',markersize=0,capsize=0,zorder=-1,
             color='0.50',ls='None',uplims=False)
-        plt.errorbar(x[uplim], 0.8*y[uplim],
+        plt.errorbar(x[uplim], y[uplim],
             yerr=[0.2*y[uplim], 0.2*y[uplim]],
+            lw=optimal_errorlinewidth,
             fmt='o',markersize=0,capsize=optimal_markersize/1.5,zorder=-1,
             color='0.50',ls='None',uplims=True)
 
@@ -498,6 +506,7 @@ def PlotSED(config,pars):
 
     #save the canvas
     #plt.grid()
+    plt.tight_layout()
     plt.savefig("%s.png" %filebase, dpi=150, facecolor='w', edgecolor='w',
             orientation='portrait', papertype=None, format=None,
             transparent=False, bbox_inch=None, pad_inches=0.1,
@@ -524,6 +533,7 @@ def PlotUL(pars,config,ULFlux,Index):
 
     #save the plot
     filebase = utils._SpecFileName(config)
+    plt.tight_layout()
     plt.savefig(filebase + '.png', dpi=150, facecolor='w', edgecolor='w',
             orientation='portrait', papertype=None, format=None,
             transparent=False, bbox_inches=None, pad_inches=0.1,