add enrico_contour

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]+' <config file name>')
+    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)

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")

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'

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 ")