Merge pull request #114 from mireianievas/master

Several changes (see below for details)

bin/enrico_download

@@ -3,7 +3,7 @@
 import os
 from optparse import OptionParser
 from enrico.data import Data
-from enrico.environ import DOWNLOAD_DIR
+from enrico.environ import DOWNLOAD_DIR, USE_FULLMISSION_SPACECRAFT
 
 parser = OptionParser(description=__doc__)
 parser.add_option("--download_data",
@@ -58,11 +58,21 @@ if options.download_data:
     from os.path import join
     weekly_dir = join(DOWNLOAD_DIR+'/weekly/photon')
     os.system('ls '+weekly_dir+'/*fits > '+join(weekly_dir,'evt.lis'))
-    liste =  open(join(weekly_dir,'evt.lis'),'r').readlines()
-    saved = open(join(DOWNLOAD_DIR,'events.lis'),"w")
+    liste    = [l.strip() for l in open(join(weekly_dir,'evt.lis'),'r').readlines()]
+    saved    = open(join(DOWNLOAD_DIR,'events.lis'),"w")
     for evt in liste:
-        saved.write(join(weekly_dir,evt)+'\n')
+        if 'p302' in evt:
+            saved.write(join(weekly_dir,evt)+'\n')
     saved.close()
+    if not USE_FULLMISSION_SPACECRAFT:
+        weekly_sc_dir = join(DOWNLOAD_DIR+'/weekly/spacecraft')
+        os.system('ls '+weekly_sc_dir+'/*fits > '+join(weekly_sc_dir,'sc.lis'))
+        liste_sc = [l.strip() for l in open(join(weekly_sc_dir,'sc.lis'),'r').readlines()]
+        saved = open(join(DOWNLOAD_DIR,'spacecraft.lis'),"w")
+        for sc in liste_sc:
+            if 'p202' in sc:
+                saved.write(join(weekly_sc_dir,sc)+'\n')
+        saved.close()
 
 if options.download_spacecraft:
     data.download(spacecraft=True, photon=False)

enrico-init.sh

@@ -15,7 +15,9 @@
 # $ init_fermi
 # $ init_enrico
 
-export ENRICO_DIR=$HOME
+if [ "$ENRICO_DIR" == "" ]; then
+    export ENRICO_DIR=$HOME
+fi
 echo "Adding Enrico to PATH and PYTHONPATH"
 export PATH=$PATH:$ENRICO_DIR/bin
 export PYTHONPATH=$ENRICO_DIR:$PYTHONPATH

enrico/RunGTlike.py

@@ -28,9 +28,40 @@ 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 is no xml file, create it and print a warning
+    if (not os.path.isfile(config['file']['xml'])):
+        mes.warning("Xml not found, creating one for the given config %s" %config['file']['xml'])
+        XmlMaker(config)
+
     Fit = SummedLikelihood.SummedLikelihood()
+    if hasKey(config,'ComponentAnalysis') == True:
+        # Create one obs instance for each component
+        configs  = [None]*4
+        Fits     = [None]*4
+        Analyses = [None]*4
+        if isKey(config['ComponentAnalysis'],'FrontBack') == 'yes':
+            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)
+                configs[k]['event']['evtype'] = fermievtypes[TYPE]
+                try:
+                    Analyses[k] = Analysis(folder, configs[k], \
+                        configgeneric=config,\
+                        tag=TYPE, verbose = verbose)
+                    if not(xmlfile ==""): Analyses[k].obs.xmlfile = xmlfile
+                    Fits[k] = Analyses[k].CreateLikeObject()
+                    Fit.addComponent(Fits[k])
+                except RuntimeError,e:
+                    if 'RuntimeError: gtltcube execution failed' in str(e):
+                        mes.warning("Event type %s is empty! Error is %s" %(TYPE,str(e)))
+            FitRunner = Analyses[0]
 
     EUnBinned = config['ComponentAnalysis']['EUnBinned']
     emintotal = float(config['energy']['emin'])

enrico/config.py

@@ -3,7 +3,7 @@
 from os.path import join
 from enrico.extern.configobj import ConfigObj, flatten_errors
 from enrico.extern.validate import Validator
-from enrico.environ import CONFIG_DIR, DOWNLOAD_DIR
+from enrico.environ import CONFIG_DIR, DOWNLOAD_DIR, USE_FULLMISSION_SPACECRAFT
 from enrico import Loggin, utils
 
 def get_config(infile, configspec=join(CONFIG_DIR, 'default.conf')):
@@ -87,7 +87,10 @@ def query_config():
 
 #    informations about the input files
     config['file'] = {}
-    config['file']['spacecraft'] = DOWNLOAD_DIR+'/lat_spacecraft_merged.fits'
+    if USE_FULLMISSION_SPACECRAFT:
+        config['file']['spacecraft'] = DOWNLOAD_DIR+'/lat_spacecraft_merged.fits'
+    else:
+        config['file']['spacecraft'] = '@'+DOWNLOAD_DIR+'/spacecraft.lis'
     ft2 = raw_input('FT2 file ['+config['file']['spacecraft']+'] : ')
     if not(ft2=='') :
       config['file']['spacecraft'] = ft2

enrico/config/default.conf

@@ -36,7 +36,7 @@ Submit = option('yes', 'no', default='no')
 	ra = float(default=0.)
 	dec = float(default=0.)
 	# supported models are PowerLaw, PowerLaw2, LogParabola, PLExpCutoff, Generic
-	spectrum = option('PowerLaw', 'PowerLaw2', 'LogParabola', 'PLExpCutoff', BrokenPowerLaw,'Generic', default='PowerLaw')
+	spectrum = option('PowerLaw', 'PowerLaw2', 'LogParabola', 'PLExpCutoff', 'BrokenPowerLaw', 'Generic', default='PowerLaw')
 	#No effect if = 0, no EBL model added to the spectral shape
 	redshift = float(default=0.)
 	ebl_model = integer(default=4, min=0, max=20)
@@ -50,6 +50,8 @@ Submit = option('yes', 'no', default='no')
 	ComputeDiffrsp = option('yes', 'no', default='yes')
 	zmax = float(default=100)
 	roicut = option('yes', 'no', default='no')
+	evtroicuts = option('yes', 'no', default='yes')
+	evttimecuts = option('yes', 'no', default='yes')
 	filter = string(default='(DATA_QUAL>0)&&(LAT_CONFIG==1)')
 
 [event]
@@ -79,13 +81,13 @@ Submit = option('yes', 'no', default='no')
 	emin = float(default=100)
 	emax = float(default=300000)
 	enumbins_per_decade = integer(default=10)	
-    decorrelation_energy = option('yes', 'no', default='no')
+	decorrelation_energy = option('yes', 'no', default='no')
 
 [fitting]
 	optimizer = option('DRMNFB', 'DRMNGB', 'MINUIT', 'NEWMINUIT', default='MINUIT')
 	ftol = float(default=1e-6)
-    # if source_ts < min_source_TS, the source is removed from the model
-    min_source_TS = float(default=1.0)
+	# if source_ts < min_source_TS, the source is removed from the model
+	min_source_TS = float(default=1.0)
 
 [model]
 	# The following options determine the xml model
@@ -165,8 +167,9 @@ Submit = option('yes', 'no', default='no')
 	NumEnergyBins = integer(default=0)
 	#Compute an UL if the TS of the sources is <TSEnergyBins
 	TSEnergyBins = float(default=9)
-    #Distribute Ebins according to the butterfly errors (optimal dist)
-    DistEbins = option('logE', 'TS', 'mix', default='logE')
+	#Distribute Ebins according to the butterfly errors (optimal dist)
+	# logE: log scaled bins, TS: similar TS, mix: in between.
+	DistEbins = option('logE', 'TS', 'mix', default='logE')
 
 
 [TSMap]

enrico/data.py

@@ -6,13 +6,13 @@
 log = logging.getLogger(__name__)
 
 # @todo: use config file for this?
-from enrico.environ import CATALOG_DIR, CATALOG, DIFFUSE_DIR, DIFFUSE_GAL, DIFFUSE_ISO_SOURCE, DIFFUSE_ISO_CLEAN
+from enrico.environ import CATALOG_DIR, CATALOG, CATALOG_8yr, DIFFUSE_DIR, DIFFUSE_GAL, DIFFUSE_ISO_SOURCE, DIFFUSE_ISO_CLEAN
 from enrico.environ import DIFFUSE_ISO_SOURCEPSF0, DIFFUSE_ISO_SOURCEPSF1, DIFFUSE_ISO_SOURCEPSF2, DIFFUSE_ISO_SOURCEPSF3
 from enrico.environ import DIFFUSE_ISO_SOURCEEDISP0, DIFFUSE_ISO_SOURCEEDISP1, DIFFUSE_ISO_SOURCEEDISP2, DIFFUSE_ISO_SOURCEEDISP3
 from enrico.environ import DIFFUSE_ISO_CLEANPSF0, DIFFUSE_ISO_CLEANPSF1, DIFFUSE_ISO_CLEANPSF2, DIFFUSE_ISO_CLEANPSF3
 from enrico.environ import DIFFUSE_ISO_CLEANEDISP0, DIFFUSE_ISO_CLEANEDISP1, DIFFUSE_ISO_CLEANEDISP2, DIFFUSE_ISO_CLEANEDISP3
 from enrico.environ import DIRS, DOWNLOAD_DIR, CATALOG_TEMPLATE_DIR, TEMPLATE_VERSION, PREPROCESSED_DIR
-from enrico.environ import WEEKLY_DIR, SPACECRAFT
+from enrico.environ import WEEKLY_DIR, SPACECRAFT, USE_FULLMISSION_SPACECRAFT
 
 #TODO: Read from default config
 default_filter = 'DATA_QUAL==1&&LAT_CONFIG==1&&ABS(ROCK_ANGLE)<52'
@@ -23,14 +23,16 @@
 HEASARC_FTP = 'ftp://heasarc.gsfc.nasa.gov/FTP/fermi/data/lat'
 WEEKLY_DIFFUSE_URL = 'ftp://heasarc.gsfc.nasa.gov/FTP/fermi/data/lat'
 CATALOG_URL = join(FSSC_URL, 'data/access/lat/4yr_catalog')
+CATALOG_URL_8yr = join(FSSC_URL, 'data/access/lat/fl8y/')
 DIFFUSE_URL = join(FSSC_URL, 'data/analysis/software/aux')
 WEEKLY_URL = join(FSSC_FTP_URL, 'weekly/photon')
+WEEKLY_SC_URL = join(FSSC_FTP_URL, 'weekly/spacecraft')
 WEEKLY_DIFFRSP_URL = join(HEASARC_FTP, 'weekly/diffuse')
-SPACECRAFT_URL = join(FSSC_FTP_URL, 'mission/spacecraft',
-                      SPACECRAFT)
+SPACECRAFT_URL = join(FSSC_FTP_URL, 'mission/spacecraft', SPACECRAFT)
 
 #        (_tag,                 _url,        _dir,        _file)
 FILES = [('CATALOG',            CATALOG_URL, CATALOG_DIR, CATALOG),
+         ('CATALOG',            CATALOG_URL_8yr, CATALOG_DIR, CATALOG_8yr),
          ('DIFFUSE_GAL',        DIFFUSE_URL, DIFFUSE_DIR, DIFFUSE_GAL),
          ('DIFFUSE_ISO_SOURCE', DIFFUSE_URL, DIFFUSE_DIR, DIFFUSE_ISO_SOURCE),
          ('DIFFUSE_ISO_SOURCEPSF0', DIFFUSE_URL, DIFFUSE_DIR, DIFFUSE_ISO_SOURCEPSF0),
@@ -125,7 +127,10 @@ def download(self, spacecraft=True, photon=True):
         os.chdir(DOWNLOAD_DIR)
         if spacecraft:
             # -m --mirror
-            cmd = 'wget -N ' + SPACECRAFT_URL
+            if USE_FULLMISSION_SPACECRAFT:
+                cmd = 'wget -N ' + SPACECRAFT_URL
+            else:
+                cmd = 'wget -m -P weekly -nH --cut-dirs=4 -np ' + WEEKLY_SC_URL
             print('EXEC: ', cmd)
             os.system(cmd)
         if photon:
@@ -232,6 +237,20 @@ def _preprocess_list(self, selection):
         files = files[:weeks]
         log.debug('Writing weeks.lis with %04d lines.' % len(files))
         open('weeks.lis', 'w').writelines(files)
+
+        if not USE_FULLMISSION_SPACECRAFT:
+            """Produce lists of weekly spacecraft list files."""
+            files = os.listdir(WEEKLY_SC_DIR)
+            # Select only fits files
+            files = [join(WEEKLY_SC_DIR, _) + '\n'
+                     for _ in files if _.endswith('.fits')]
+            # Sort chronologically
+            files.sort()
+            # Select only a subset of weeks
+            weeks = self.SELECTIONS[selection]
+            files = files[:weeks]
+            log.debug('Writing weeks.lis with %04d lines.' % len(files))
+            open('weeks_sc.lis', 'w').writelines(files)
 
     def _preprocess_gtselect(self, evclass, emin):
         """Run gtselect"""
@@ -258,7 +277,10 @@ def _preprocess_gtmktime(self):
         """Run gtmktime"""
         from gt_apps import maketime as tool
         self._set_common_tool_options(tool)
-        tool['scfile'] = join(DOWNLOAD_DIR, SPACECRAFT)
+        if USE_FULLMISSION_SPACECRAFT:
+            tool['scfile'] = join(DOWNLOAD_DIR, SPACECRAFT)
+        else:
+            tool['scfile'] = '@weeks_sc.lis'
         tool['sctable'] = 'SC_DATA'
         tool['filter'] = default_filter
         tool['roicut'] = 'no'
@@ -278,7 +300,10 @@ def _preprocess_gtltcube(self):
         self._set_common_tool_options(tool)
         tool['evfile'] = 'gtmktime.fits'
         tool['evtable'] = 'EVENTS'
-        tool['scfile'] = join(DOWNLOAD_DIR, SPACECRAFT)
+        if USE_FULLMISSION_SPACECRAFT:
+            tool['scfile'] = join(DOWNLOAD_DIR, SPACECRAFT)
+        else:
+            tool['scfile'] = '@weeks_sc.lis'
         tool['sctable'] = 'SC_DATA'
         tool['outfile'] = 'gtltcube.fits'
         tool['dcostheta'] = 0.025

enrico/energybin.py

@@ -4,7 +4,7 @@
 from enrico.constants import EbinPath
 from enrico.submit import call
 from enrico.config import get_config
-from enrico import utils,Loggin
+from enrico import utils, Loggin
 
 def ChangeModel(comp, E1, E2, name, Pref, Gamma):
     """Change the spectral model of a source called name
@@ -115,7 +115,7 @@ def PrepareEbin(Fit, FitRunner,sedresult=None):
         xmltag_list = ["_FRONT","_BACK"]
         mes.info("Splitting Front/Back events")
     elif config['ComponentAnalysis']['PSF'] == "yes":
-        xmltag_list = ["_PSF0","_PSF1","_PSF2"]
+        xmltag_list = ["_PSF0","_PSF1","_PSF2","_PSF3"]
         mes.info("Splitting PSF events")
     elif config['ComponentAnalysis']['EDISP'] == "yes":
         xmltag_list = ["_EDISP0","_EDISP1","_EDISP2","_EDISP3"]
@@ -126,7 +126,7 @@ def PrepareEbin(Fit, FitRunner,sedresult=None):
         E = utils.GetE0(ener[ibin + 1],ener[ibin])
         mes.info("Submitting # "+str(ibin)+" at energy "+str(E))
         #Update the model for the bin
-        for  comp,xmltag in zip(Fit.components, xmltag_list):
+        for comp,xmltag in zip(Fit.components, xmltag_list):
             NewFitObject = ChangeModel(comp, ener[ibin], ener[ibin + 1], srcname, Pref[ibin] ,Gamma[ibin])
             Xmlname = (config['out'] + "/" + srcname + "_" + str(ibin) +xmltag+ ".xml")
 

enrico/environ.py

@@ -25,11 +25,13 @@
 DIFFUSE_DIR = os.environ.get('FERMI_DIFFUSE_DIR', '')
 DOWNLOAD_DIR = os.environ.get('FERMI_DOWNLOAD_DIR', '')
 WEEKLY_DIR = ''
+WEEKLY_SC_DIR = ''
 if DOWNLOAD_DIR :
   WEEKLY_DIR = join(DOWNLOAD_DIR, 'weekly/photon')
+  WEEKLY_SC_DIR = join(DOWNLOAD_DIR, 'weekly/spacecraft')
 PREPROCESSED_DIR = os.environ.get('FERMI_PREPROCESSED_DIR', '')
 CONFIG_DIR = join(os.path.dirname(__file__), 'config')
-
+USE_FULLMISSION_SPACECRAFT = bool(os.environ.get('USE_FULLMISSION_SPACECRAFT','False')=='True')
 
 try :
     from enrico.extern.odict import OrderedDict
@@ -40,13 +42,15 @@
                    PREPROCESSED_DIR=PREPROCESSED_DIR,
                    DOWNLOAD_DIR=DOWNLOAD_DIR,
                    WEEKLY_DIR=WEEKLY_DIR,
+                   WEEKLY_SC_DIR=WEEKLY_SC_DIR,
                    CONFIG_DIR=CONFIG_DIR,
                    ENRICO_DIR=ENRICO_DIR)
 except :
     DIRS = {}
 
 # File names
 CATALOG_VERSION = '16'
+CATALOG_8yr_VERSION = '5'
 TEMPLATE_VERSION = '15'
 if CATALOG_DIR :
   CATALOG_TEMPLATE_DIR = join(CATALOG_DIR, 'Extended_archive_v%s/Templates'% TEMPLATE_VERSION)
@@ -56,6 +60,7 @@
     pass
 
 CATALOG = 'gll_psc_v%s.fit' % CATALOG_VERSION
+CATALOG_8yr = 'gll_psc_8year_v%s.fit' % CATALOG_8yr_VERSION
 DIFFUSE_GAL = 'gll_iem_v06.fits'
 DIFFUSE_ISO_SOURCE = 'iso_P8R2_SOURCE_V6_v06.txt'
 DIFFUSE_ISO_SOURCE = 'iso_P8R2_SOURCE_V6_v06.txt'

enrico/fitmaker.py

@@ -104,7 +104,8 @@ def CreateLikeObject(self):
                               irfs=self.obs.irfs)
             Fit = UnbinnedAnalysis(Obs, self.obs.xmlfile,
                                    optimizer=self.config['fitting']['optimizer'])
-
+
+	# Fix this, EBL absorbed models use LogParabola with b=0 instead of PowerLaw, we may want to allow fixed shape for that case
         if float(self.config['Spectrum']['FrozenSpectralIndex']>0) and self.config['target']['spectrum'] == "PowerLaw":
             parameters = dict()
             parameters['Index']  = -float(self.config['Spectrum']['FrozenSpectralIndex'])