diff --git a/bin/enrico_contour b/bin/enrico_contour new file mode 100755 index 0000000..eea808c --- /dev/null +++ b/bin/enrico_contour @@ -0,0 +1,30 @@ +#!/usr/bin/env python +"""Make a lightcurve""" +import sys +import numpy as np +import enrico.scan as scan +from enrico.config import get_config +from enrico import Loggin +mes = Loggin.Message() +try: + infile = sys.argv[1] +except: + print('Usage: '+sys.argv[0]+' ') + mes.error('Config file not found.') + +if len(sys.argv)==2 : +#read an config file alone. If not working, try to read an ascii file with different conf file + try : + liste = np.genfromtxt(sys.argv[1],dtype="str",unpack=True) + for inf in liste: + mes.info("work on the config file "+inf) + config = get_config(inf) + scan.Scan(config) + except : + config = get_config(infile) + scan.Contour(config) +else: + for inf in sys.argv[1:]: + mes.info("work on the config file "+inf) + config = get_config(inf) + scan.Scan(config) diff --git a/enrico/config/default.conf b/enrico/config/default.conf index 921d6b7..885c61f 100644 --- a/enrico/config/default.conf +++ b/enrico/config/default.conf @@ -182,3 +182,7 @@ Submit = option('yes', 'no', default='no') rad = float(default=1) # list of sources srclist = string(default="") + +[Contours] + parname1 = string(default="Prefactor") + parname2 = string(default="Index") diff --git a/enrico/environ.py b/enrico/environ.py index 6854b14..6cfb962 100755 --- a/enrico/environ.py +++ b/enrico/environ.py @@ -47,6 +47,10 @@ TEMPLATE_VERSION = '15' if CATALOG_DIR : CATALOG_TEMPLATE_DIR = join(CATALOG_DIR, 'Extended_archive_v%s/Templates'% TEMPLATE_VERSION) + try : + os.mkdir(join(CATALOG_DIR, 'Extended_archive_v%s'% TEMPLATE_VERSION)) + except: + pass CATALOG = 'gll_psc_v%s.fit' % CATALOG_VERSION DIFFUSE_GAL = 'gll_iem_v05_rev1.fit' diff --git a/enrico/scan.py b/enrico/scan.py index 31f5b7e..b478472 100644 --- a/enrico/scan.py +++ b/enrico/scan.py @@ -5,18 +5,22 @@ import enrico.constants as cst import RunGTlike import ROOT -import numpy,os,string +import numpy,os,string,array +from enrico import Loggin -def MakeScan(Fit,spectrum,par,bmin,bmax,opt): - N=30 + +def MakeScan(Fit,spectrum,par,bmin,bmax,opt,N=100): Param = numpy.zeros(N) loglike = numpy.zeros(N) for i in xrange(N): Param[i] = bmin + (bmax-bmin)*i/(N-1.) spectrum.getParam(par).setValue(Param[i]) - spectrum.getParam(par).setFree(0) - loglike[i] = Fit.fit(0,covar=False,optimizer=opt) - +# spectrum.getParam(par).setFree(0) +# loglike[i] = Fit.fit(0,covar=False,optimizer=opt) +# Fit.optimize(0) + Fit.logLike.syncParams() + loglike[i] = -Fit.logLike.value() + return Param,loglike def Scan(config): @@ -60,3 +64,97 @@ def Scan(config): cres.Print(config["out"]+"/"+cst.ScanPath+ "/Scan_"+par+".C") cres.Print(config["out"]+"/"+cst.ScanPath+ "/Scan_"+par+".png") +def Contour(config): + ROOT.gROOT.SetBatch(ROOT.kTRUE) +# cres = ROOT.TCanvas("Contour") + config["Spectrum"]["FitsGeneration"] = "no" + parname1 = config["Contours"]["parname1"] + parname2 = config["Contours"]["parname2"] + + FitRunner,Fit = RunGTlike.GenAnalysisObjects(config) + spectrum = Fit[FitRunner.obs.srcname].funcs['Spectrum'] + + ParName = spectrum.paramNames + + mes = Loggin.Message() + mes.info("Computing Contours for "+parname1+" and "+parname2) + + + ### Check part !!!! + findpar2 = findpar1 = False + for par in ParName : #Loop over the parameters to check + if par == parname1: + findpar1 = True + if not(spectrum.getParam(par).isFree()): + mes.error(parname1+" is not a free parameter") + if par == parname2: + findpar2 = True + if not(spectrum.getParam(par).isFree()): + mes.error(parname2+" is not a free parameter") + + if not(findpar1): + mes.error(parname1+" is not a valid parameter") + if not(findpar2): + mes.error(parname2+" is not a valid parameter") + + bestloglike = Fit.fit(0,covar=False,optimizer=config['fitting']['optimizer']) + print spectrum + print "Min LogLikelihood =",bestloglike + + ## get values + ParValue1 = spectrum.getParam(parname1).value() + ParError1 = spectrum.getParam(parname1).error() + bmin1,bmax1 = spectrum.getParam(parname1).getBounds() + + bmin1 = max(bmin1,ParValue1-20*ParError1) + bmax1 = min(bmax1,ParValue1+20*ParError1) + + ParValue2 = spectrum.getParam(parname2).value() + ParError2 = spectrum.getParam(parname2).error() + bmin2,bmax2 = spectrum.getParam(parname2).getBounds() + + bmin2 = max(bmin2,ParValue2-20*ParError2) + bmax2 = min(bmax2,ParValue2+20*ParError2) + + N = 100 + param2 = numpy.zeros(N) + loglike = ROOT.TH2F("loglike","Contours (68%, 95%, 99%)",N,bmin1,bmax1,N,bmin2,bmax2) + spectrum.getParam(parname2).setFree(0) + + mes.info("Boundaries for "+parname1+" ["+str(bmin1)+","+str(bmax1)+"]") + mes.info("Boundaries for "+parname2+" ["+str(bmin2)+","+str(bmax2)+"]") + + for i in xrange(N): + + param2[i] = bmin2 + (bmax2-bmin2)*i/(N-1.) + + spectrum.getParam(parname2).setValue(param2[i]) + + param1,ll = MakeScan(Fit,spectrum,parname1,bmin1,bmax1,config['fitting']['optimizer'],N) + + for j in xrange(N): + loglike.Fill(param1[j],param2[i],ll[j]) + + os.system("mkdir -p "+config["out"]+"/"+cst.ScanPath) + cres = ROOT.TCanvas("Contours") + loglike.SetMinimum(bestloglike); + loglike.SetMaximum(bestloglike+3); + loglike.SetXTitle(parname1); + loglike.SetYTitle(parname2); + + loglike.SetStats(000) + loglike.SetContour(3) + loglike.SetContourLevel(0,bestloglike+0.5) + loglike.SetContourLevel(1,bestloglike+4./2.) + loglike.SetContourLevel(2,bestloglike+6.63/2.) + loglike.Draw("CONT1"); + + tgrres = ROOT.TGraphErrors(2,array.array('f',[ParValue1,ParValue1]),array.array('f',[ParValue2,ParValue2]),array.array('f',[ParError1,ParError1]),array.array('f',[ParError2,ParError2])) + tgrres.Draw(".pz") + cres.Print(config["out"]+"/"+cst.ScanPath+ "/Contours_"+parname1+"_"+parname2+".eps") + cres.Print(config["out"]+"/"+cst.ScanPath+ "/Contours_"+parname1+"_"+parname2+".C") + cres.Print(config["out"]+"/"+cst.ScanPath+ "/Contours_"+parname1+"_"+parname2+".png") + + + + mes.success("Scan ")