Skip to content

Commit

Permalink
some improvement
Browse files Browse the repository at this point in the history
  • Loading branch information
@davidsanchez
davidsanchez committed Feb 14, 2018
1 parent 302dc0c commit 946a4bf
Show file tree
Hide file tree
Showing 11 changed files with 265 additions and 232 deletions.
1 change: 1 addition & 0 deletions enrico-init.sh
View file
Expand Up @@ -15,6 +15,7 @@
# $ init_fermi
# $ init_enrico

export ENRICO_DIR=/Users/davidsanchez/Soft/enrico
echo "Adding Enrico to PATH and PYTHONPATH"
export PATH=$PATH:$ENRICO_DIR/bin
export PYTHONPATH=$ENRICO_DIR:$PYTHONPATH
Expand Down
206 changes: 70 additions & 136 deletions enrico/RunGTlike.py
View file
@@ -1,11 +1,14 @@
#!/usr/bin/env python
import os,os.path,math
from enrico import utils
from enrico.gtfunction import Observation
from enrico.fitmaker import FitMaker
import Loggin
import SummedLikelihood
from enrico.extern.configobj import ConfigObj
from utils import hasKey, isKey, typeirfs

def Analysis(folder, config, configgeneric=None, tag="", convtyp='-1', verbose = 1):
import os
from enrico import utils
from enrico.gtfunction import Observation
from enrico.fitmaker import FitMaker
import Loggin

mes = Loggin.Message()
""" run an analysis"""
Expand All @@ -14,154 +17,85 @@ def Analysis(folder, config, configgeneric=None, tag="", convtyp='-1', verbose =
utils._log('SUMMARY: ' + tag)
Obs.printSum()


FitRunner = FitMaker(Obs, config)##Class
if config['Spectrum']['FitsGeneration'] == 'yes':
FitRunner.FirstSelection(configgeneric) #Generates fits files for the coarse selection
FitRunner.GenerateFits() #Generates fits files for the rest of the products
return FitRunner

def GenAnalysisObjects(config, verbose = 1, xmlfile =""):
import os
import os.path
import math
import SummedLikelihood
from enrico.xml_model import XmlMaker
from enrico.extern.configobj import ConfigObj
from utils import hasKey, isKey
import Loggin

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]

elif isKey(config['ComponentAnalysis'],'PSF') == 'yes':
from enrico.data import fermievtypes
mes.info("Breaking the analysis in PSF 0,1,2,3.")
# Clone the configs
# Set Summed Likelihood to True
config['Spectrum']['SummedLike'] = 'yes'
for k,TYPE in enumerate(["PSF0", "PSF1", "PSF2", "PSF3"]):
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]

elif isKey(config['ComponentAnalysis'],'EDISP') == 'yes':
from enrico.data import fermievtypes
mes.info("Breaking the analysis in EDISP 0,1,2,3.")
# Clone the configs
# Set Summed Likelihood to True
config['Spectrum']['SummedLike'] = 'yes'
for k,TYPE in enumerate(["EDISP0", "EDISP1", "EDISP2", "EDISP3"]):
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]

elif isKey(config['ComponentAnalysis'],'EUnBinned')>=0:

EUnBinned = config['ComponentAnalysis']['EUnBinned']
emintotal = float(config['energy']['emin'])
emaxtotal = float(config['energy']['emax'])

evtnum = [config["event"]["evtype"]] #for std analysis
evtold = evtnum[0] #for std analysis

# Create one obs instance for each component
if isKey(config['ComponentAnalysis'],'FrontBack') == 'yes':
evtnum = [1, 2]
if isKey(config['ComponentAnalysis'],'PSF') == 'yes':
evtnum = [4,8,16,32]
if isKey(config['ComponentAnalysis'],'EDISP') == 'yes':
evtnum = [64,128,256,521]
oldxml = config['file']['xml']
for k,evt in enumerate(evtnum):
config['event']['evtype'] = evt
config["file"]["xml"] = oldxml.replace(".xml","_"+typeirfs[evt]+".xml").replace("_.xml",".xml")

if EUnBinned>emintotal and EUnBinned<emaxtotal:
mes.info("Breaking the analysis in Binned (low energy) and Unbinned (high energies)")
# Set Summed Likelihood to True
config['Spectrum']['SummedLike'] = 'yes'
EUnBinned = config['ComponentAnalysis']['EUnBinned']
emintotal = float(config['energy']['emin'])
emaxtotal = float(config['energy']['emax'])
# Run the following analysis depending on the case
# this is general, no matter if we are in the total fit or the Ebin #N fit.
if EUnBinned<=emintotal: analysestorun = ["highE"]
elif EUnBinned>=emaxtotal: analysestorun = ["lowE"]
else: analysestorun = ["lowE","highE"]
oldxml = config['file']['xml']
analysestorun = ["lowE","highE"]

for k,TYPE in enumerate(analysestorun):
configs[k] = ConfigObj(config)
tag = TYPE
if typeirfs[evt] != "" : tag += "_"+typeirfs[evt]# handle name of fits file

# Tune parameters
if TYPE is "lowE":
configs[k]['energy']['emin'] = emintotal
configs[k]['energy']['emax'] = min(config['energy']['emax'],EUnBinned)
configs[k]['analysis']['likelihood'] = "binned"
configs[k]['analysis']['ComputeDiffrsp'] = "no"
config['energy']['emin'] = emintotal
config['energy']['emax'] = min(config['energy']['emax'],EUnBinned)
config['analysis']['likelihood'] = "binned"
config['analysis']['ComputeDiffrsp'] = "no"
elif TYPE is "highE":
configs[k]['energy']['emin'] = max(config['energy']['emin'],EUnBinned)
configs[k]['energy']['emax'] = emaxtotal
configs[k]['analysis']['likelihood'] = "unbinned"
configs[k]['analysis']['ComputeDiffrsp'] = "yes"
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:
raise
if 'RuntimeError: gtltcube execution failed' in str(e):
mes.warning("Event type %s is empty! Error is %s" %(TYPE,str(e)))
FitRunner = Analyses[0]
config['energy']['emin'] = max(config['energy']['emin'],EUnBinned)
config['energy']['emax'] = emaxtotal
config['analysis']['likelihood'] = "unbinned"
config['analysis']['ComputeDiffrsp'] = "yes"

Analyse = Analysis(folder, config, \
configgeneric=config,\
tag=TYPE,\
verbose=verbose)


Fit_component = Analyse.CreateLikeObject()
Fit.addComponent(Fit_component)
FitRunner = Analyse
FitRunner.obs.Emin = emintotal
FitRunner.obs.Emax = emaxtotal

try:
FitRunner.config
except NameError:
# Create one obs instance
FitRunner = Analysis(folder, config, tag="", verbose = verbose)
Fit.addComponent(FitRunner.CreateLikeObject())
if not(xmlfile ==""):
FitRunner.obs.xmlfile = xmlfile
else:
Analyse = Analysis(folder, config, \
configgeneric=config,\
tag=typeirfs[evt], verbose = verbose)

# if not(xmlfile ==""): Analyse.obs.xmlfile = xmlfile
Fit_component = Analyse.CreateLikeObject()
Fit.addComponent(Fit_component)
FitRunner = Analyse

config["event"]["evtype"] = evtold
FitRunner.config = config

return FitRunner,Fit

def run(infile):
Expand All @@ -179,6 +113,7 @@ def run(infile):
FitRunner,Fit = GenAnalysisObjects(config)
# create all the fit files and run gtlike
FitRunner.PerformFit(Fit)
sedresult = None

#plot the SED and model map if possible and asked
if float(config['UpperLimit']['TSlimit']) < Fit.Ts(config['target']['name']):
Expand Down Expand Up @@ -215,16 +150,15 @@ def run(infile):

if config['Spectrum']['ResultPlots'] == 'yes' :
outXml = utils._dump_xml(config)
if config['Spectrum']['SummedLike'] != 'yes':
# the possiblity of making the model map is checked inside the function
FitRunner.ModelMap(outXml)


# the possibility of making the model map is checked inside the function
FitRunner.ModelMap(outXml)

# Make energy bins by running a *new* analysis
Nbin = config['Ebin']['NumEnergyBins']
energybin.RunEbin(folder,Nbin,Fit,FitRunner)

energybin.RunEbin(folder,Nbin,Fit,FitRunner,sedresult)

del(sedresult)
del(Result)
del(FitRunner)

Expand Down
13 changes: 7 additions & 6 deletions enrico/config/default.conf
View file
Expand Up @@ -36,13 +36,13 @@ 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', '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)
# 0=Kneiske, 1=Primack05, 2=Kneiske_HighUV, 3=Stecker05, 4=Franceschini, 5=Finke, 6=Gilmore, 10=Dominguez11
# full list http://www-glast.stanford.edu/cgi-bin/viewcvs/celestialSources/eblAtten/eblAtten/EblAtten.h?view=markup
fit_tau = option('yes', 'no', default='yes')
fit_tau = option('yes', 'no', default='no')

[analysis]
# General analysis options
Expand Down Expand Up @@ -79,7 +79,7 @@ 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='yes')
decorrelation_energy = option('yes', 'no', default='no')

[fitting]
optimizer = option('DRMNFB', 'DRMNGB', 'MINUIT', 'NEWMINUIT', default='MINUIT')
Expand All @@ -106,8 +106,6 @@ Submit = option('yes', 'no', default='no')
ResultPlots = option('yes', 'no', default='yes')
#Freeze the spectral index of the source. Has no implication if 0 (Left free to vary)
FrozenSpectralIndex = float(default=0, min=0, max=5)
#Use the summed likelihood method
SummedLike = option('yes', 'no', default='no')


[UpperLimit]
Expand Down Expand Up @@ -167,6 +165,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')


[TSMap]
#Re-fit before computing the TS map
Expand Down Expand Up @@ -204,7 +205,7 @@ Submit = option('yes', 'no', default='no')
# binned and unbinned analysis should be split
# disabled if -1
EUnBinned = float(default=-1)
# Divide the analysis in PSF classes
# Divide the analysis in FrontBack classes
FrontBack = option('yes', 'no', default='no')
# Divide the analysis in PSF classes
PSF = option('yes', 'no', default='no')
Expand Down
24 changes: 12 additions & 12 deletions enrico/data.py
View file
Expand Up @@ -52,18 +52,18 @@
('DIFFUSE_ISO_CLEAN', DIFFUSE_URL, DIFFUSE_DIR, DIFFUSE_ISO_CLEAN)]


fermievtypes = {\
'FRONT': 1,
'BACK': 2,
'PSF0': 4,
'PSF1': 8,
'PSF2': 16,
'PSF3': 32,
'EDISP0': 64,
'EDISP1': 128,
'EDISP2': 256,
'EDISP3': 512,
}
# fermievtypes = {\
# 'FRONT': 1,
# 'BACK': 2,
# 'PSF0': 4,
# 'PSF1': 8,
# 'PSF2': 16,
# 'PSF3': 32,
# 'EDISP0': 64,
# 'EDISP1': 128,
# 'EDISP2': 256,
# 'EDISP3': 512,
# }


def check_dirs():
Expand Down

0 comments on commit 946a4bf

Please sign in to comment.