ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/run.py
Revision: 1.2
Committed: Tue May 8 10:26:31 2012 UTC (13 years ago) by peller
Content type: text/x-python
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +0 -0 lines
State: FILE REMOVED
Log Message: 
test

File Contents

# User Rev Content
1 peller 1.1 #!/usr/bin/env python
2     import sys
3     import os
4     import ROOT
5     import shutil
6     from ROOT import TFile
7     from array import array
8     from init import *
9     #from initMuMu import *
10     from math import sqrt
11     import random
12     from copy import copy
13     #suppres the EvalInstace conversion warning bug
14     import warnings
15     warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
16     from ConfigParser import SafeConfigParser
17     from samplesinfo import sample
18     from mvainfos import mvainfo
19     import pickle
20     from progbar import progbar
21    
22    
23     for i in range(0,len(backgroundFiles)): backgroundFiles[i] = prefix + backgroundFiles[i]
24     for i in range(0,len(signalFiles)): signalFiles[i] = prefix + signalFiles[i]
25     for i in range(0,len(dataFiles)): dataFiles[i] = prefix + dataFiles[i]
26     #for i in range(0,len(InFiles0)): InFiles0[i] = Preprefix + InFiles0[i]
27     #for i in range(0,len(InFiles1)): InFiles1[i] = Preprefix + InFiles1[i]
28     #for i in range(0,len(InFiles2)): InFiles2[i] = Preprefix + InFiles2[i]
29    
30     #add ehm together
31     jobs = copy(backgroundFiles)
32     jobs.append(signalFiles[0])
33     legenden = backname
34     legenden.append(signame[0])
35     xsecs = xsec
36     xsecs.append(signal_xsec)
37     if sys.argv[1] == 'stack' or 'pie': treeCut = treeCutPlot
38     else: treeCut = treeCutMVA
39    
40    
41    
42    
43     #CONFIGURE
44    
45     #load config
46     config = SafeConfigParser()
47     config.read('./config')
48    
49     #get locations:
50     Wdir=config.get('Directories','Wdir')
51    
52    
53     #systematics
54     systematics=config.get('systematics','systematics')
55     systematics=systematics.split(' ')
56    
57     #TreeVar Array
58     MVA_Vars={}
59     for systematic in systematics:
60     MVA_Vars[systematic]=config.get('treeVars',systematic)
61     MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
62    
63    
64    
65    
66    
67    
68    
69    
70    
71    
72    
73    
74    
75     ##############################
76     # #
77     # VHbb Analysis Stuff #
78     # ETHZ, Philipp Eller #
79     # #
80     ##############################
81    
82    
83     #*********************HELPERS*******************************
84     def getTree(path,job):
85     Tree = ROOT.TChain(treeName)
86     Tree.Add("%s/%s.root" %(path,job))
87     Tree.SetDirectory(0)
88     #nEntries = Tree.GetEntries()
89     return Tree
90    
91     #NEW
92     def getTree2(job):
93     Tree = ROOT.TChain(job.tree)
94     Tree.Add(job.getpath())
95     Tree.SetDirectory(0)
96     return Tree
97    
98     def getScale2(job):
99     input = TFile.Open(job.getpath())
100     CountWithPU = input.Get("CountWithPU")
101     CountWithPU2011B = input.Get("CountWithPU2011B")
102     #print lumi*xsecs[i]/hist.GetBinContent(1)
103     return float(job.lumi)*float(job.xsec)/(0.46502*CountWithPU.GetBinContent(1)+0.53498*CountWithPU2011B.GetBinContent(1))*1/float(job.split)
104    
105     def getScale(path,job):
106     input = TFile.Open("%s/%s.root" %(path,job))
107     CountWithPU = input.Get("CountWithPU")
108     CountWithPU2011B = input.Get("CountWithPU2011B")
109     i=jobs.index(job)
110     return lumi*xsecs[i]/(0.46502*CountWithPU.GetBinContent(1)+0.53498*CountWithPU2011B.GetBinContent(1))
111    
112    
113     #NEW
114     def AddSystematics(path):
115    
116     infofile = open(path+'/samples.info','r')
117     info = pickle.load(infofile)
118     infofile.close()
119     os.mkdir(path+'/sys')
120    
121     for job in info:
122    
123     if job.type != 'DATA':
124     print '\t - %s' %(job.name)
125    
126     input = TFile.Open(job.getpath(),'read')
127     Count = input.Get("Count")
128     CountWithPU = input.Get("CountWithPU")
129     CountWithPU2011B = input.Get("CountWithPU2011B")
130     tree = input.Get(job.tree)
131     nEntries = tree.GetEntries()
132    
133     job.addpath('/sys')
134     output = ROOT.TFile(job.getpath(), 'RECREATE')
135     newtree = tree.CloneTree(0)
136     job.SYS = ['Nominal','JER_up','JER_down','JES_up','JES_down','beff_up','beff_down','bmis_up','bmis_down']
137    
138    
139     hJ0 = ROOT.TLorentzVector()
140     hJ1 = ROOT.TLorentzVector()
141    
142     #JER branches
143     hJet_pt_JER_up = array('f',[0]*2)
144     newtree.Branch('hJet_pt_JER_up',hJet_pt_JER_up,'hJet_pt_JER_up[2]/F')
145     hJet_pt_JER_down = array('f',[0]*2)
146     newtree.Branch('hJet_pt_JER_down',hJet_pt_JER_down,'hJet_pt_JER_down[2]/F')
147     hJet_e_JER_up = array('f',[0]*2)
148     newtree.Branch('hJet_e_JER_up',hJet_e_JER_up,'hJet_e_JER_up[2]/F')
149     hJet_e_JER_down = array('f',[0]*2)
150     newtree.Branch('hJet_e_JER_down',hJet_e_JER_down,'hJet_e_JER_down[2]/F')
151     H_JER = array('f',[0]*4)
152     newtree.Branch('H_JER',H_JER,'mass_up:mass_down:pt_up:pt_down/F')
153    
154     #JES branches
155     hJet_pt_JES_up = array('f',[0]*2)
156     newtree.Branch('hJet_pt_JES_up',hJet_pt_JES_up,'hJet_pt_JES_up[2]/F')
157     hJet_pt_JES_down = array('f',[0]*2)
158     newtree.Branch('hJet_pt_JES_down',hJet_pt_JES_down,'hJet_pt_JES_down[2]/F')
159     hJet_e_JES_up = array('f',[0]*2)
160     newtree.Branch('hJet_e_JES_up',hJet_e_JES_up,'hJet_e_JES_up[2]/F')
161     hJet_e_JES_down = array('f',[0]*2)
162     newtree.Branch('hJet_e_JES_down',hJet_e_JES_down,'hJet_e_JES_down[2]/F')
163     H_JES = array('f',[0]*4)
164     newtree.Branch('H_JES',H_JES,'mass_up:mass_down:pt_up:pt_down/F')
165    
166     #Add training Flag
167     EventForTraining = array('f',[0])
168     newtree.Branch('EventForTraining',EventForTraining,'EventForTraining/F')
169    
170     iter=0
171    
172     for entry in range(0,nEntries):
173     tree.GetEntry(entry)
174    
175     #fill training flag
176     iter+=1
177     if (iter%2==0):
178     EventForTraining=1
179     else:
180     EventForTraining=0
181    
182     hJet_pt0 = tree.hJet_pt[0]
183     hJet_pt1 = tree.hJet_pt[1]
184     hJet_eta0 = tree.hJet_eta[0]
185     hJet_eta1 = tree.hJet_eta[1]
186     hJet_genPt0 = tree.hJet_genPt[0]
187     hJet_genPt1 = tree.hJet_genPt[1]
188     hJet_e0 = tree.hJet_e[0]
189     hJet_e1 = tree.hJet_e[1]
190     hJet_phi0 = tree.hJet_phi[0]
191     hJet_phi1 = tree.hJet_phi[1]
192     hJet_JECUnc0 = tree.hJet_JECUnc[0]
193     hJet_JECUnc1 = tree.hJet_JECUnc[1]
194    
195    
196     for updown in ['up','down']:
197     #JER
198     if updown == 'up':
199     inner = 0.06
200     outer = 0.1
201     if updown == 'down':
202     inner = -0.06
203     outer = -0.1
204     #Calculate
205     if abs(hJet_eta0)<1.1: res0 = inner
206     else: res0 = outer
207     if abs(hJet_eta1)<1.1: res1 = inner
208     else: res1 = outer
209     rPt0 = hJet_pt0 + (hJet_pt0-hJet_genPt0)*res0
210     rPt1 = hJet_pt1 + (hJet_pt1-hJet_genPt1)*res1
211     rE0 = hJet_e0*rPt0/hJet_pt0
212     rE1 = hJet_e1*rPt1/hJet_pt1
213     hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
214     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
215     #Set
216     if updown == 'up':
217     hJet_pt_JER_up[0]=rPt0
218     hJet_pt_JER_up[1]=rPt1
219     hJet_e_JER_up[0]=rE0
220     hJet_e_JER_up[1]=rE1
221     H_JER[0]=(hJ0+hJ1).M()
222     H_JER[2]=(hJ0+hJ1).Pt()
223     if updown == 'down':
224     hJet_pt_JER_down[0]=rPt0
225     hJet_pt_JER_down[1]=rPt1
226     hJet_e_JER_down[0]=rE0
227     hJet_e_JER_down[1]=rE1
228     H_JER[1]=(hJ0+hJ1).M()
229     H_JER[3]=(hJ0+hJ1).Pt()
230    
231     #JES
232     if updown == 'up':
233     variation=1
234     if updown == 'down':
235     variation=-1
236     #calculate
237     rPt0 = hJet_pt0*(1+variation*hJet_JECUnc0)
238     rPt1 = hJet_pt1*(1+variation*hJet_JECUnc1)
239     rE0 = hJet_e0*(1+variation*hJet_JECUnc0)
240     rE1 = hJet_e1*(1+variation*hJet_JECUnc1)
241     hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
242     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
243     #Fill
244     if updown == 'up':
245     hJet_pt_JES_up[0]=rPt0
246     hJet_pt_JES_up[1]=rPt1
247     hJet_e_JES_up[0]=rE0
248     hJet_e_JES_up[1]=rE1
249     H_JES[0]=(hJ0+hJ1).M()
250     H_JES[2]=(hJ0+hJ1).Pt()
251     if updown == 'down':
252     hJet_pt_JES_down[0]=rPt0
253     hJet_pt_JES_down[1]=rPt1
254     hJet_e_JES_down[0]=rE0
255     hJet_e_JES_down[1]=rE1
256     H_JES[1]=(hJ0+hJ1).M()
257     H_JES[3]=(hJ0+hJ1).Pt()
258    
259     newtree.Fill()
260    
261     newtree.Write()
262     Count.Write()
263     CountWithPU.Write()
264     CountWithPU2011B.Write()
265     output.Close()
266    
267     else: #(is data)
268    
269     shutil.copy(job.getpath(),path+'/sys')
270     job.addpath('/sys')
271    
272     infofile = open(path+'/sys'+'/samples.info','w')
273     pickle.dump(info,infofile)
274     infofile.close()
275    
276     #NEW
277     def Addcut(path,addpath,Samplecut,Datacut):
278    
279     infofile = open(path+'/samples.info','r')
280     info = pickle.load(infofile)
281     infofile.close()
282     os.mkdir(path+addpath)
283     addprefix=''
284    
285     Samplecut=config.get('Cuts',Samplecut)
286     Datacut=config.get('Cuts',Datacut)
287    
288    
289     for job in info:
290     if job.type != 'DATA':
291     cut = Samplecut
292     print '\t - %s' %(job.name)
293     copytree2(job,addpath,addprefix,cut)
294     job.addpath(addpath)
295     job.addtreecut(cut)
296     job.addprefix(addprefix)
297     job.addcomment('added cut ' + cut)
298    
299     if job.type == 'DATA':
300     cut = Datacut
301     print '\t - %s' %(job.name)
302     copytree2(job,addpath,addprefix,cut)
303     job.addpath(addpath)
304     job.addtreecut(cut)
305     job.addprefix(addprefix)
306     job.addcomment('added cut ' + cut)
307    
308     infofile = open(path+addpath+'/samples.info','w')
309     pickle.dump(info,infofile)
310     infofile.close()
311    
312     def Addsinglecut(path,name,prefix,cut):
313    
314     infofile = open(path+'/samples.info','r')
315     info = pickle.load(infofile)
316     infofile.close()
317    
318     for job in info:
319     if job.name == name:
320     print '\t - %s' %(job.name)
321     copytree2(job,'',prefix,cut)
322     job.addtreecut(cut)
323     job.addcomment('added cut ' + cut)
324    
325     infofile = open(path+'/samples.info','w')
326     pickle.dump(info,infofile)
327     infofile.close()
328    
329     def AddFile(path,name,newname,prefix,cut):
330    
331     infofile = open(path+'/samples.info','r')
332     info = pickle.load(infofile)
333     infofile.close()
334    
335     for job in info:
336     if job.name == name:
337     print '\t - %s' %(job.name)
338     copytree2(job,'',prefix,cut)
339     job2 = copy(job)
340     job2.addtreecut(cut)
341     #job2.addprefix(prefix)
342     job2.addcomment('added cut ' + cut)
343     job2.name=newname
344    
345     info.append(job2)
346     infofile = open(path+'/samples.info','w')
347     pickle.dump(info,infofile)
348     infofile.close()
349    
350     def getHistoFromTree(path,job,scale):
351     Tree = getTree(path,job)
352     hTree = ROOT.TH1F(job,job,nBins,xMin,xMax)
353     if scale !=0:
354     Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,job,nBins,xMin,xMax),'%s*(%s)' %(weightF,treeCut), "goff")
355     else:
356     Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,job,nBins,xMin,xMax),'(%s)' %(treeCut), "goff")
357     hTree = ROOT.gDirectory.Get(job)
358     hTree.SetDirectory(0)
359     if scale !=0:
360     ScaleFactor = getScale(treePath,job)
361     if ScaleFactor != 0:
362     hTree.Scale(ScaleFactor)
363     return hTree
364    
365     def getHistoFromTree2(job,options):
366     Tree = getTree2(job)
367     treeVar=options[0]
368     name=job.name
369     title=job.plotname()
370     nBins=int(options[3])
371     xMin=float(options[4])
372     xMax=float(options[5])
373     if job.type != 'DATA':
374     treeCut=config.get('Cuts',options[7])
375     elif job.type == 'DATA':
376     treeCut=config.get('Cuts',options[8])
377     weightF=config.get('Weights','weightF')
378     hTree = ROOT.TH1F('%s'%name,'%s'%title,nBins,xMin,xMax)
379     hTree.Sumw2()
380     if job.type != 'DATA':
381     Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'%s*(%s)' %(weightF,treeCut), "goff,e")
382     else:
383     Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'(%s)' %(treeCut), "goff,e")
384     hTree = ROOT.gDirectory.Get(name)
385     hTree.SetDirectory(0)
386     #print job.name + ' Sumw2', hTree.GetEntries()
387    
388     if job.type != 'DATA':
389     ScaleFactor = getScale2(job)
390     if ScaleFactor != 0:
391     hTree.Scale(ScaleFactor)
392    
393     hTree.SetFillStyle(1001)
394     hTree.SetFillColor(job.plotcolor())
395     hTree.SetLineWidth(1)
396     print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
397    
398    
399     return hTree, job.plotname()
400    
401     def copytreeetc(path,job,Nprefix,Acut):
402     input = TFile.Open("%s/%s%s.root" %(skim_path,Preprefix,job))
403     Count = input.Get("Count")
404     CountWithPU = input.Get("CountWithPU")
405     CountWithPU2011B = input.Get("CountWithPU2011B")
406     inputTree = input.Get(treeName)
407     nEntries = inputTree.GetEntries()
408     output = TFile.Open("%s/%s%s%s.root" %(path,Preprefix,Nprefix,job),'recreate')
409     print 'copy file: ' + job
410     outputTree = inputTree.CopyTree(Acut)
411     kEntries = outputTree.GetEntries()
412     #factor = kEntries/nEntries
413     print "\t before cuts\t %s" %nEntries
414     print "\t survived\t %s" %kEntries
415     #print factor
416     #print "\t Factor for Scaling is %s" %factor
417     outputTree.AutoSave()
418     #Count.Scale(factor)
419     Count.Write()
420     CountWithPU.Write()
421     #CountWithPU.Scale(factor)
422     CountWithPU2011B.Write()
423     #CountWithPU2011B.Scale(factor)
424     input.Close()
425     output.Close()
426    
427     #NEW
428     def copytree2(job,addpath,addprefix,addcut):
429     input = TFile.Open(job.getpath(),'read')
430     Count = input.Get("Count")
431     CountWithPU = input.Get("CountWithPU")
432     CountWithPU2011B = input.Get("CountWithPU2011B")
433     inputTree = input.Get(job.tree)
434     nEntries = inputTree.GetEntries()
435     output = TFile.Open("%s%s/%s%s%s.root" %(job.path,addpath,addprefix,job.prefix,job.identifier()),'recreate')
436     print '\t\tcopy file: ' + job.name + ' with cut ' + addcut + ' to ' + addpath
437     outputTree = inputTree.CopyTree(addcut)
438     kEntries = outputTree.GetEntries()
439     print "\t before cuts\t %s" %nEntries
440     print "\t survived\t %s" %kEntries
441     outputTree.AutoSave()
442     Count.Write()
443     CountWithPU.Write()
444     CountWithPU2011B.Write()
445     input.Close()
446     output.Close()
447    
448    
449     def CutCopy(mode):
450    
451     if mode == 'all':
452     path = treePath
453     for job in InFiles0:
454     copytreeetc(path,job,prefix0[0],cut0[0])#+'&'+treeCutMVA)
455     for i in range(2):
456     copytreeetc(path,InFiles1[0],prefix1[i],cut1[i])#+'&'+treeCutMVA)
457     copytreeetc(path,InFiles2[0],prefix2[i],cut2[i])#+'&'+treeCutMVA)
458    
459     if mode == 'Zbb':
460     path = treePath + '/Zbb'
461     for job in InFiles0:
462     copytreeetc(path,job,prefix0[0],cut0[0]+'&'+treeCutZbb)
463     for i in range(2):
464     copytreeetc(path,InFiles1[0],prefix1[i],cut1[i]+'&'+treeCutZbb)
465     copytreeetc(path,InFiles2[0],prefix2[i],cut2[i]+'&'+treeCutZbb)
466     for job in dataFiles0:
467     copytreeetc(path,job,prefix0[0],cut0[0]+'&'+treeCutZbb)
468     if mode == 'Zlight':
469     path = treePath + '/Zlight'
470     for job in InFiles0:
471     copytreeetc(path,job,prefix0[0],cut0[0]+'&'+treeCutZlight)
472     for i in range(2):
473     copytreeetc(path,InFiles1[0],prefix1[i],cut1[i]+'&'+treeCutZlight)
474     copytreeetc(path,InFiles2[0],prefix2[i],cut2[i]+'&'+treeCutZlight)
475     for job in dataFiles0:
476     copytreeetc(path,job,prefix0[0],cut0[0]+'&'+treeCutZlight)
477     if mode == 'Top':
478     path = treePath + '/Top'
479     for job in InFiles0:
480     copytreeetc(path,job,prefix0[0],cut0[0]+'&'+treeCutTop)
481     for i in range(2):
482     copytreeetc(path,InFiles1[0],prefix1[i],cut1[i]+'&'+treeCutTop)
483     copytreeetc(path,InFiles2[0],prefix2[i],cut2[i]+'&'+treeCutTop)
484     for job in dataFiles0:
485     copytreeetc(path,job,prefix0[0],cut0[0]+'&'+treeCutTop)
486     if mode == 'Signal':
487     path = treePath + '/Signal'
488     for job in InFiles0:
489     copytreeetc(path,job,prefix0[0],cut0[0]+'&'+treeCutSignal)
490     for i in range(2):
491     copytreeetc(path,InFiles1[0],prefix1[i],cut1[i]+'&'+treeCutSignal)
492     copytreeetc(path,InFiles2[0],prefix2[i],cut2[i]+'&'+treeCutSignal)
493     for job in dataFiles0:
494     copytreeetc(path,job,prefix0[0],cut0[0]+'&'+treeCutSignal)
495    
496     if mode == 'data':
497     path = treePath #+ '/test'
498     for job in dataFiles0:
499     copytreeetc(path,job,prefix0[0],cut0[0]+'&'+triggerFlags[0])#+'&'+treeCutMVA)
500    
501     if mode == 'test':
502     path = treePath +'/test'
503     for job in InFiles0:
504     copytreeetc(path,job,prefix0[0],cut0[0]+'&'+cut_even+'&'+treeCutMVA)
505     for i in range(2):
506     copytreeetc(path,InFiles1[0],prefix1[i],cut1[i]+'&'+cut_even+'&'+treeCutMVA)
507     copytreeetc(path,InFiles2[0],prefix2[i],cut2[i]+'&'+cut_even+'&'+treeCutMVA)
508    
509     if mode == 'train':
510     path = treePath +'/train'
511     for job in InFiles0:
512     copytreeetc(path,job,prefix0[0],cut0[0]+'&'+cut_odd+'&'+treeCutMVA)
513     for i in range(2):
514     copytreeetc(path,InFiles1[0],prefix1[i],cut1[i]+'&'+cut_odd+'&'+treeCutMVA)
515     copytreeetc(path,InFiles2[0],prefix2[i],cut2[i]+'&'+cut_odd+'&'+treeCutMVA)
516    
517     def createComparison():
518     ROOT.gROOT.SetStyle("Plain")
519     c = ROOT.TCanvas(title,title, 800, 600)
520     c.SetLogy()
521     allStack = ROOT.THStack(title,title)
522     histos = []
523     for job in jobs:
524     hTemp = getHistoFromTree(treePath,job)
525     histos.append(hTemp)
526     l = ROOT.TLegend(0.55, 0.82, 0.93, 0.93)
527     for i in range(0,len(histos)):
528     histos[i].SetFillStyle(0)
529     histos[i].SetLineColor(i+2)
530     histos[i].SetLineWidth(2)
531     allStack.Add(histos[i])
532     l.AddEntry(histos[i],legenden[i],'l')
533     allStack.Draw("HISTNOSTACK")
534     allStack.GetXaxis().SetTitle(xTitle)
535     allStack.GetYaxis().SetTitle(yTitle)
536     allStack.GetXaxis().SetRangeUser(xMin,xMax)
537     l.Draw()
538     name = '%s/Comparison/%s.png' %(plotPath,title)
539     c.Print(name)
540    
541     def plot():
542     print 'ok, i make plots for you now...'
543     histos = []
544     for job in jobs:
545     hTemp = getHistoFromTree(treePath,job,1)
546     histos.append(hTemp)
547     datas = []
548     for job in dataFiles:
549     hTemp = getHistoFromTree(treePath,job,0)
550     datas.append(hTemp)
551     createStack(histos,datas,plotPath,title,treeVar,xMin,xMax,nBins,xTitle,yTitle)
552    
553     def allplots():
554     print 'ok, i make all plots for you now...'
555     for i in range(0,len(treeVars)-1):
556     global treeVar
557     treeVar = treeVars[i]
558     xTitle = treeVar
559     title = set + '_' + treeVar
560     global xMin
561     xMin = xMinS[i]
562     global xMax
563     xMax = xMaxS[i]
564     global nBins
565     nBins = nBinsS[i]
566     histos = []
567     for job in jobs:
568     hTemp = getHistoFromTree(treePath,job,1)
569     histos.append(hTemp)
570     datas = []
571     for job in dataFiles:
572     hTemp = getHistoFromTree(treePath,job,0)
573     datas.append(hTemp)
574     createStack(histos,datas,plotPath,title,treeVar,xMin,xMax,nBins,xTitle,yTitle)
575    
576     def createStack(histos,datas,plotPath,title,treeVar,xMin,xMax,nBins,xTitle,yTitle):
577     #*********************STACK*******************************
578     print '*******************'
579     print 'now i am working on ' + title
580     ROOT.gROOT.SetStyle("Plain")
581     c = ROOT.TCanvas(title,title, 800, 600)
582     allStack = ROOT.THStack(title,title)
583     l = ROOT.TLegend(0.68, 0.63, 0.88, 0.88)
584     MC_integral=0
585     for i in range(0,len(histos)):
586     histos[i].SetFillStyle(1001)
587     histos[i].SetFillColor(color[i])
588     histos[i].SetLineWidth(1)
589     #histos[i].SetLineColor(color[i])
590     print "\t%s integral:" %legenden[i]
591     print histos[i].Integral()
592     MC_integral=MC_integral+histos[i].Integral()
593     print "MC integral:"
594     print MC_integral
595     histos[0].Add(histos[1],1)
596     del histos[1]
597     histos[1].Add(histos[2],1)
598     del histos[2]
599     for i in range(2):
600     histos[2].Add(histos[3],1)
601     del histos[3]
602     for i in range(5):
603     histos[4].Add(histos[5],1)
604     del histos[5]
605     k=len(histos)
606     for j in range(0,k):
607     #print histos[j].GetBinContent(1)
608     i=k-j-1
609     allStack.Add(histos[i])
610     l.AddEntry(histos[j],legende[j],'F')
611     d1 = ROOT.TH1F('d1','d1',nBins,xMin,xMax)
612    
613     for i in range(0,len(datas)):
614     d1.Add(datas[i],1)
615     print "data integral:"
616     print d1.Integral()
617     print d1.GetEntries()
618     l.AddEntry(d1,datalegend,'PL')
619     allStack.SetTitle()
620     allStack.Draw("")
621     allStack.GetXaxis().SetTitle(xTitle)
622     allStack.GetYaxis().SetTitle(yTitle)
623     allStack.GetXaxis().SetRangeUser(xMin,xMax)
624     allStack.GetYaxis().SetRangeUser(0,20000)
625     Ymax = max(allStack.GetMaximum(),d1.GetMaximum())*1.2
626     allStack.SetMaximum(Ymax)
627     allStack.SetMinimum(0.01)
628     c.Update()
629     ROOT.gPad.SetLogy()
630     ROOT.gPad.SetTicks(1,1)
631     allStack.Draw("")
632     d1.SetMarkerStyle(21)
633     d1.Draw("P0,E1,X0,same")
634     l.SetFillColor(0)
635     l.SetBorderSize(0)
636     l.Draw()
637     name = '%s/Stack/%s.png' %(plotPath,title)
638     c.Print(name)
639    
640    
641     def treeStack(path,var,data):
642     #*********************STACK*******************************
643    
644     plot=config.get('Plot',var)
645    
646     infofile = open(path+'/samples.info','r')
647     info = pickle.load(infofile)
648     infofile.close()
649    
650     options = plot.split(',')
651     name=options[1]
652     title = options[2]
653     nBins=int(options[3])
654     xMin=float(options[4])
655     xMax=float(options[5])
656    
657     setup=config.get('Plot','setup')
658     setup=setup.split(',')
659    
660    
661     print '\nProducing Plot of %s\n'%title
662    
663    
664     histos = []
665     typs = []
666     datas = []
667     datatyps =[]
668    
669     for job in info:
670     if job.name != 'DATA':
671     hTemp, typ = getHistoFromTree2(job,options)
672     histos.append(hTemp)
673     typs.append(typ)
674     elif job.name in data:
675     hTemp, typ = getHistoFromTree2(job,options)
676     datas.append(hTemp)
677     datatyps.append(typ)
678    
679    
680    
681    
682     ROOT.gROOT.SetStyle("Plain")
683     c = ROOT.TCanvas(name,title, 800, 600)
684     allStack = ROOT.THStack(name,title)
685     l = ROOT.TLegend(0.68, 0.63, 0.88, 0.88)
686     MC_integral=0
687     MC_entries=0
688    
689     for histo in histos:
690     MC_integral+=histo.Integral()
691     #MC_entries+=histo.GetEntries()
692     print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
693     #flow = MC_entries-MC_integral
694     #if flow > 0:
695     # print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
696    
697     #ORDER AND ADD TOGETHER
698    
699     ordnung=[]
700     ordnungtyp=[]
701     num=[0]*len(setup)
702     for i in range(0,len(setup)):
703     for j in range(0,len(histos)):
704     if typs[j] == setup[i]:
705     num[i]+=1
706     ordnung.append(histos[j])
707     ordnungtyp.append(typs[j])
708    
709     del histos
710     del typs
711    
712     histos=ordnung
713     typs=ordnungtyp
714    
715     for k in range(0,len(num)):
716     for m in range(0,num[k]):
717     if m > 0:
718     histos[k].Add(histos[k+1],1)
719     del histos[k+1]
720     del typs[k+1]
721    
722     k=len(histos)
723     for j in range(0,k):
724     #print histos[j].GetBinContent(1)
725     i=k-j-1
726     allStack.Add(histos[i])
727     l.AddEntry(histos[j],typs[j],'F')
728    
729    
730     d1 = ROOT.TH1F('d1','d1',nBins,xMin,xMax)
731    
732     datatitle=''
733     for i in range(0,len(datas)):
734     d1.Add(datas[i],1)
735     datatitle=datatitle+ ' + '+datatyps[i]
736     print "\033[1;32m\n\tDATA integral = %s\033[1;m"%d1.Integral()
737     flow = d1.GetEntries()-d1.Integral()
738     if flow > 0:
739     print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
740     l.AddEntry(d1,datatitle,'PL')
741     allStack.SetTitle()
742     allStack.Draw("")
743     allStack.GetXaxis().SetTitle(title)
744     allStack.GetYaxis().SetTitle('Counts')
745     allStack.GetXaxis().SetRangeUser(xMin,xMax)
746     allStack.GetYaxis().SetRangeUser(0,20000)
747     Ymax = max(allStack.GetMaximum(),d1.GetMaximum())*1.2
748     allStack.SetMaximum(Ymax)
749     allStack.SetMinimum(0.1)
750     c.Update()
751     if config.get('Plot','logy') == '1':
752     ROOT.gPad.SetLogy()
753     ROOT.gPad.SetTicks(1,1)
754     allStack.Draw("")
755     d1.SetMarkerStyle(21)
756     d1.Draw("P0,E1,X0,same")
757     l.SetFillColor(0)
758     l.SetBorderSize(0)
759     l.Draw()
760    
761     name = '%s/%s' %(config.get('Directories','plotpath'),options[6])
762     c.Print(name)
763    
764     def treeCompare(path,var):
765     #*********************STACK*******************************
766    
767    
768     infofile = open(path+'/samples.info','r')
769     info = pickle.load(infofile)
770     infofile.close()
771    
772     plot=config.get('Compare',var)
773     options = plot.split(',')
774     name=options[1]
775     title = options[2]
776     nBins=int(options[3])
777     xMin=float(options[4])
778     xMax=float(options[5])
779    
780    
781     bkgs=config.get('Compare','BKG')
782     bkgs=bkgs.split(' ')
783    
784     sigs=config.get('Compare','SIG')
785     sigs=sigs.split(' ')
786    
787    
788     bkgsA = []
789     bkgsB = []
790     sigsA = []
791     sigsB = []
792    
793     optionsA=copy(options)
794     optionsB=copy(options)
795    
796     optionsA[7]=config.get('Compare','cutA')
797     optionsB[7]=config.get('Compare','cutB')
798    
799    
800    
801    
802     for job in info:
803     if job.name in bkgs:
804     hTemp, typ = getHistoFromTree2(job,optionsA)
805     bkgsA.append(hTemp)
806     hTemp, typ = getHistoFromTree2(job,optionsB)
807     bkgsB.append(hTemp)
808     if job.name in sigs:
809     hTemp, typ = getHistoFromTree2(job,optionsA)
810     sigsA.append(hTemp)
811     hTemp, typ = getHistoFromTree2(job,optionsB)
812     sigsB.append(hTemp)
813    
814    
815     ROOT.gROOT.SetStyle("Plain")
816     c = ROOT.TCanvas(name,title, 800, 600)
817     l = ROOT.TLegend(0.68, 0.7, 0.88, 0.88)
818    
819     ROOT.gStyle.SetOptStat(0)
820     ROOT.gROOT.ForceStyle
821    
822     for i in range(1,len(sigsA)):
823     sigsA[0].Add(sigsA[i])
824     sigsB[0].Add(sigsB[i])
825    
826    
827     for i in range(1,len(bkgsA)):
828     bkgsA[0].Add(bkgsA[i])
829     bkgsB[0].Add(bkgsB[i])
830    
831     ScaleFactor=1/sigsA[0].Integral()
832     sigsA[0].Scale(ScaleFactor)
833     ScaleFactor=1/sigsB[0].Integral()
834     sigsB[0].Scale(ScaleFactor)
835     ScaleFactor=1/bkgsA[0].Integral()
836     bkgsA[0].Scale(ScaleFactor)
837     ScaleFactor=1/bkgsB[0].Integral()
838     bkgsB[0].Scale(ScaleFactor)
839    
840    
841     sigsA[0].SetLineColor(2)
842     sigsB[0].SetLineColor(1)
843     bkgsA[0].SetLineColor(4)
844     bkgsB[0].SetLineColor(1)
845    
846     sigsA[0].SetFillColor(0)
847     sigsB[0].SetMarkerColor(2)
848     bkgsA[0].SetFillColor(0)
849     bkgsB[0].SetMarkerColor(4)
850    
851    
852     sigsA[0].SetLineWidth(2)
853     sigsB[0].SetLineWidth(1)
854     bkgsA[0].SetLineWidth(2)
855     bkgsB[0].SetLineWidth(1)
856    
857    
858     sigsA[0].SetFillStyle(3000)
859     sigsB[0].SetFillStyle(3345)
860     bkgsA[0].SetFillStyle(3000)
861     bkgsB[0].SetFillStyle(3354)
862    
863     l.AddEntry(sigsA[0],'SIG %s'%optionsA[7],'L')
864     l.AddEntry(sigsB[0],'SIG %s'%optionsB[7],'PL')
865     l.AddEntry(bkgsA[0],'BKG %s'%optionsA[7],'L')
866     l.AddEntry(bkgsB[0],'BKG %s'%optionsB[7],'PL')
867    
868    
869     maximum=max(sigsA[0].GetMaximum(),sigsB[0].GetMaximum(),bkgsA[0].GetMaximum(),bkgsB[0].GetMaximum())
870    
871    
872     sigsA[0].SetTitle("Comparison EE/MM")
873     sigsA[0].Draw("hist")
874     sigsA[0].GetXaxis().SetTitle(title)
875     sigsA[0].GetYaxis().SetTitle('Normalized Scale')
876     sigsA[0].GetXaxis().SetRangeUser(xMin,xMax)
877     sigsA[0].GetYaxis().SetRangeUser(0,maximum*1.1)
878     #c.Update()
879     #if config.get('Plot','logy') == '1':
880     # ROOT.gPad.SetLogy()
881     ROOT.gPad.SetTicks(1,1)
882    
883    
884     l.SetFillColor(0)
885     l.SetBorderSize(0)
886     l.Draw()
887     print sigsA[0].GetEntries()
888     print sigsB[0].GetEntries()
889     print bkgsA[0].GetEntries()
890     print bkgsB[0].GetEntries()
891    
892    
893     sigsA[0].Draw("hist,same")
894     sigsB[0].SetMarkerStyle(21)
895     #sigsB[0].Sumw2()
896     sigsB[0].Draw("P0,same")
897     bkgsA[0].Draw("hist,same")
898     bkgsB[0].SetMarkerStyle(21)
899     #bkgsB[0].Sumw2()
900     bkgsB[0].Draw("P0,same")
901    
902    
903    
904    
905     name = '%s/%s' %(config.get('Directories','plotpath'),options[6])
906     c.Print(name)
907    
908     #NEW TRAINING
909     def newTraining(run,gui):
910    
911     #CONFIG
912     #factory
913     factoryname=config.get('factory','factoryname')
914     factorysettings=config.get('factory','factorysettings')
915     #MVA
916     MVAtype=config.get(run,'MVAtype')
917     MVAname=config.get(run,'MVAname')
918     MVAsettings=config.get(run,'MVAsettings')
919     fnameOutput = Wdir +'/weights/'+factoryname+'_'+MVAname+'.root'
920     #locations
921     Tpath=config.get(run,'Tpath')
922     Epath=config.get(run,'Epath')
923    
924     #signals
925     signals=config.get(run,'signals')
926     signals=signals.split(' ')
927     #backgrounds
928     backgrounds=config.get(run,'backgrounds')
929     backgrounds=backgrounds.split(' ')
930    
931     treeVarSet=config.get(run,'treeVarSet')
932    
933     #variables
934     #TreeVar Array
935     MVA_Vars={}
936     MVA_Vars['Nominal']=config.get(treeVarSet,'Nominal')
937     MVA_Vars['Nominal']=MVA_Vars['Nominal'].split(' ')
938     #Spectators:
939     spectators=config.get(treeVarSet,'spectators')
940     spectators=spectators.split(' ')
941    
942     #TRAINING samples
943     infofile = open(Tpath+'/samples.info','r')
944     Tinfo = pickle.load(infofile)
945     infofile.close()
946    
947     #Evaluate samples
948     infofile = open(Epath+'/samples.info','r')
949     Einfo = pickle.load(infofile)
950     infofile.close()
951    
952     #Workdir
953     workdir=ROOT.gDirectory.GetPath()
954    
955    
956     #load TRAIN trees
957     Tbackgrounds = []
958     TbScales = []
959     Tsignals = []
960     TsScales = []
961    
962     for job in Tinfo:
963     if job.name in signals:
964     Tsignal = getTree2(job)
965     ROOT.gDirectory.Cd(workdir)
966     TsScale = getScale2(job)
967     Tsignals.append(Tsignal)
968     TsScales.append(TsScale)
969    
970     if job.name in backgrounds:
971     Tbackground = getTree2(job)
972     ROOT.gDirectory.Cd(workdir)
973     TbScale = getScale2(job)
974     Tbackgrounds.append(Tbackground)
975     TbScales.append(TbScale)
976    
977     #load EVALUATE trees
978     Ebackgrounds = []
979     EbScales = []
980     Esignals = []
981     EsScales = []
982    
983     for job in Einfo:
984     if job.name in signals:
985     Esignal = getTree2(job)
986     ROOT.gDirectory.Cd(workdir)
987     EsScale = getScale2(job)
988     Esignals.append(Esignal)
989     EsScales.append(EsScale)
990    
991     if job.name in backgrounds:
992     Ebackground = getTree2(job)
993     ROOT.gDirectory.Cd(workdir)
994     EbScale = getScale2(job)
995     Ebackgrounds.append(Ebackground)
996     EbScales.append(EbScale)
997    
998     output = ROOT.TFile.Open(fnameOutput, "RECREATE")
999     factory = ROOT.TMVA.Factory(factoryname, output, factorysettings)
1000    
1001     #set input trees
1002     for i in range(len(Tsignals)):
1003    
1004     factory.AddSignalTree(Tsignals[i], TsScales[i], ROOT.TMVA.Types.kTraining)
1005     factory.AddSignalTree(Esignals[i], EsScales[i], ROOT.TMVA.Types.kTesting)
1006    
1007     for i in range(len(Tbackgrounds)):
1008     if (Tbackgrounds[i].GetEntries()>0):
1009     factory.AddBackgroundTree(Tbackgrounds[i], TbScales[i], ROOT.TMVA.Types.kTraining)
1010    
1011     if (Ebackgrounds[i].GetEntries()>0):
1012     factory.AddBackgroundTree(Ebackgrounds[i], EbScales[i], ROOT.TMVA.Types.kTesting)
1013    
1014    
1015     for var in MVA_Vars['Nominal']:
1016     factory.AddVariable(var,'D') # add the variables
1017     for var in spectators:
1018     factory.AddSpectator(var,'D') #add specators
1019    
1020     #Execute TMVA
1021     factory.SetSignalWeightExpression(weightF)
1022     factory.Verbose()
1023     factory.BookMethod(MVAtype,MVAname,MVAsettings)
1024     factory.TrainAllMethods()
1025     factory.TestAllMethods()
1026     factory.EvaluateAllMethods()
1027     output.Write()
1028    
1029     #WRITE INFOFILE
1030     infofile = open(Wdir+'/weights/'+factoryname+'_'+MVAname+'.info','w')
1031     info=mvainfo(MVAname)
1032     info.factoryname=factoryname
1033     info.factorysettings=factorysettings
1034     info.MVAtype=MVAtype
1035     info.MVAsettings=MVAsettings
1036     info.weightfilepath=Wdir+'/weights'
1037     info.Tpath=Tpath
1038     info.Epath=Epath
1039     info.varset=treeVarSet
1040     info.vars=MVA_Vars['Nominal']
1041     info.spectators=spectators
1042     pickle.dump(info,infofile)
1043     infofile.close()
1044    
1045     # open the TMVA Gui
1046     if gui == 'gui':
1047     ROOT.gROOT.ProcessLine( ".L TMVAGui.C")
1048     ROOT.gROOT.ProcessLine( "TMVAGui(\"%s\")" % fnameOutput )
1049     ROOT.gApplication.Run()
1050    
1051    
1052     #*********************EVALUATE*******************************
1053     def evaluate(run,Apath): #fnameOutput='training.root', split=0.5
1054    
1055     #CONFIG
1056     #factory
1057     factoryname=config.get('factory','factoryname')
1058     #MVA
1059     MVAname=config.get(run,'MVAname')
1060     #print Wdir+'/weights/'+factoryname+'_'+MVAname+'.info'
1061     MVAinfofile = open(Wdir+'/weights/'+factoryname+'_'+MVAname+'.info','r')
1062     MVAinfo = pickle.load(MVAinfofile)
1063     treeVarSet=MVAinfo.varset
1064     #variables
1065     #TreeVar Array
1066     MVA_Vars={}
1067     for systematic in systematics:
1068     MVA_Vars[systematic]=config.get(treeVarSet,systematic)
1069     MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
1070     #Spectators:
1071     spectators=config.get(treeVarSet,'spectators')
1072     spectators=spectators.split(' ')
1073     #progbar quatsch
1074     longe=40
1075     #Workdir
1076     workdir=ROOT.gDirectory.GetPath()
1077    
1078     os.mkdir(Apath+'/'+run)
1079    
1080    
1081     #Book TMVA readers: MVAlist=["MMCC_bla","CC5050_bla"]
1082     reader=ROOT.TMVA.Reader("!Color:!Silent" )
1083    
1084     #define variables and specatators
1085     MVA_var_buffer = []
1086     for i in range(len( MVA_Vars['Nominal'])):
1087     MVA_var_buffer.append(array( 'f', [ 0 ] ))
1088     reader.AddVariable( MVA_Vars['Nominal'][i],MVA_var_buffer[i])
1089     MVA_spectator_buffer = []
1090     for i in range(len(spectators)):
1091     MVA_spectator_buffer.append(array( 'f', [ 0 ] ))
1092     reader.AddSpectator(spectators[i],MVA_spectator_buffer[i])
1093     #Load raeder
1094     reader.BookMVA(MVAinfo.MVAname,MVAinfo.getweightfile())
1095     #--> Now the MVA is booked
1096    
1097     #Apply samples
1098     infofile = open(Apath+'/samples.info','r')
1099     Ainfo = pickle.load(infofile)
1100     infofile.close()
1101    
1102     #eval
1103     for job in Ainfo:
1104     job.addcomment('Added MVA %s'%MVAinfo.MVAname)
1105     #MCs
1106    
1107    
1108    
1109    
1110     input = TFile.Open(job.getpath(),'read')
1111     Count = input.Get("Count")
1112     CountWithPU = input.Get("CountWithPU")
1113     CountWithPU2011B = input.Get("CountWithPU2011B")
1114     tree = input.Get(job.tree)
1115     nEntries = tree.GetEntries()
1116    
1117    
1118    
1119    
1120     #tree = getTree2(job)
1121     ROOT.gDirectory.Cd(workdir)
1122     #nEntries=tree.GetEntries()
1123    
1124    
1125    
1126    
1127     if job.type != 'DATA':
1128    
1129    
1130     MVA_formulas={}
1131     for systematic in systematics:
1132     #print '\t\t - ' + systematic
1133     MVA_formulas[systematic]=[]
1134     #create TTreeFormulas
1135     for j in range(len( MVA_Vars['Nominal'])):
1136     MVA_formulas[systematic].append(ROOT.TTreeFormula("MVA_formula%s_%s"%(j,systematic),MVA_Vars[systematic][j],tree))
1137     job.addpath('/%s'%run)
1138     outfile = ROOT.TFile(job.getpath(), 'RECREATE')
1139     newtree = tree.CloneTree(0)
1140     #Setup Branches
1141     MVA = array('f',[0]*9)
1142     newtree.Branch(MVAinfo.MVAname,MVA,'nominal:JER_up:JER_down:JES_up:JES_down:beff_up:beff_down:bmis_up:bmis_down/F')
1143     print '\n--> ' + job.name +':'
1144     #progbar setup
1145     if nEntries >= longe:
1146     step=int(nEntries/longe)
1147     long=longe
1148     else:
1149     long=nEntries
1150     step = 1
1151     bar=progbar(long)
1152    
1153     #Fill event by event:
1154     for entry in range(0,nEntries):
1155     if entry % step == 0:
1156     bar.move()
1157     #load entry
1158     tree.GetEntry(entry)
1159     for systematic in systematics:
1160     for j in range(len( MVA_Vars['Nominal'])):
1161     MVA_var_buffer[j][0] = MVA_formulas[systematic][j].EvalInstance()
1162     MVA[systematics.index(systematic)] = reader.EvaluateMVA(MVAinfo.MVAname)
1163     #Fill:
1164     newtree.Fill()
1165    
1166     newtree.Write()
1167     newtree.Write()
1168     Count.Write()
1169     CountWithPU.Write()
1170     CountWithPU2011B.Write()
1171     outfile.Close()
1172    
1173     #DATA
1174     if job.type == 'DATA':
1175    
1176     #MVA Formulas
1177     MVA_formulas_Nominal = []
1178     #create TTreeFormulas
1179     for j in range(len( MVA_Vars['Nominal'])):
1180     MVA_formulas_Nominal.append(ROOT.TTreeFormula("MVA_formula%s_Nominal"%j, MVA_Vars['Nominal'][j],tree))
1181     job.addpath('/%s'%run)
1182     outfile = ROOT.TFile(job.getpath(), 'RECREATE')
1183     newtree = tree.CloneTree(0)
1184     #Setup Branches
1185     MVA = array('f',[0])
1186     newtree.Branch(MVAinfo.MVAname,MVA,'nominal/F')
1187     #progbar
1188     print '\n--> ' + job.name +':'
1189     if nEntries >= longe:
1190     step=int(nEntries/longe)
1191     long=longe
1192     else:
1193     long=nEntries
1194     step = 1
1195     bar=progbar(long)
1196    
1197     #Fill event by event:
1198     for entry in range(0,nEntries):
1199     if entry % step == 0:
1200     bar.move()
1201     #load entry
1202     tree.GetEntry(entry)
1203     #nominal:
1204     for j in range(len( MVA_Vars['Nominal'])):
1205     MVA_var_buffer[j][0] = MVA_formulas_Nominal[j].EvalInstance()
1206     MVA[0]= discr = reader.EvaluateMVA(MVAinfo.MVAname)
1207     newtree.Fill()
1208     newtree.Write()
1209     newtree.Write()
1210     Count.Write()
1211     CountWithPU.Write()
1212     CountWithPU2011B.Write()
1213     outfile.Close()
1214    
1215     print '\n'
1216     infofile = open(Apath+'/'+run+'/samples.info','w')
1217     pickle.dump(Ainfo,infofile)
1218     infofile.close()
1219    
1220     ######################
1221     #Evaluate multi: Must Have same treeVars!!!
1222     def evalMulti(Apath,arglist): #fnameOutput='training.root', split=0.5
1223     MVAlist=arglist.split(',')
1224    
1225     #CONFIG
1226     #factory
1227     factoryname=config.get('factory','factoryname')
1228     #MVA
1229     MVAnames=[]
1230     for MVA in MVAlist:
1231     print MVA
1232     MVAnames.append(config.get(MVA,'MVAname'))
1233     #print Wdir+'/weights/'+factoryname+'_'+MVAname+'.info'
1234     #MVAinfofiles=[]
1235     MVAinfos=[]
1236     for MVAname in MVAnames:
1237     MVAinfofile = open(Wdir+'/weights/'+factoryname+'_'+MVAname+'.info','r')
1238     MVAinfos.append(pickle.load(MVAinfofile))
1239     MVAinfofile.close()
1240    
1241     treeVarSet=MVAinfos[0].varset
1242     #variables
1243     #TreeVar Array
1244     MVA_Vars={}
1245     for systematic in systematics:
1246     MVA_Vars[systematic]=config.get(treeVarSet,systematic)
1247     MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
1248     #Spectators:
1249     spectators=config.get(treeVarSet,'spectators')
1250     spectators=spectators.split(' ')
1251     #progbar quatsch
1252     longe=40
1253     #Workdir
1254     workdir=ROOT.gDirectory.GetPath()
1255     os.mkdir(Apath+'/MVAout')
1256    
1257     #Book TMVA readers: MVAlist=["MMCC_bla","CC5050_bla"]
1258     readers=[]
1259     for MVA in MVAlist:
1260     readers.append(ROOT.TMVA.Reader("!Color:!Silent"))
1261    
1262     #define variables and specatators
1263     MVA_var_buffer = []
1264     for i in range(len( MVA_Vars['Nominal'])):
1265     MVA_var_buffer.append(array( 'f', [ 0 ] ))
1266     for reader in readers:
1267     reader.AddVariable( MVA_Vars['Nominal'][i],MVA_var_buffer[i])
1268     MVA_spectator_buffer = []
1269     for i in range(len(spectators)):
1270     MVA_spectator_buffer.append(array( 'f', [ 0 ] ))
1271     for reader in readers:
1272     reader.AddSpectator(spectators[i],MVA_spectator_buffer[i])
1273     #Load raeder
1274     for i in range(0,len(readers)):
1275     readers[i].BookMVA(MVAinfos[i].MVAname,MVAinfos[i].getweightfile())
1276     #--> Now the MVA is booked
1277    
1278     #Apply samples
1279     infofile = open(Apath+'/samples.info','r')
1280     Ainfo = pickle.load(infofile)
1281     infofile.close()
1282    
1283     #eval
1284     for job in Ainfo:
1285     for MVAinfo in MVAinfos:
1286     job.addcomment('Added MVA %s'%MVAinfo.MVAname)
1287     #MCs
1288    
1289     input = TFile.Open(job.getpath(),'read')
1290     Count = input.Get("Count")
1291     CountWithPU = input.Get("CountWithPU")
1292     CountWithPU2011B = input.Get("CountWithPU2011B")
1293     tree = input.Get(job.tree)
1294     nEntries = tree.GetEntries()
1295    
1296     #tree = getTree2(job)
1297     ROOT.gDirectory.Cd(workdir)
1298     #nEntries=tree.GetEntries()
1299    
1300     if job.type != 'DATA':
1301    
1302     MVA_formulas={}
1303     for systematic in systematics:
1304     #print '\t\t - ' + systematic
1305     MVA_formulas[systematic]=[]
1306     #create TTreeFormulas
1307     for j in range(len( MVA_Vars['Nominal'])):
1308     MVA_formulas[systematic].append(ROOT.TTreeFormula("MVA_formula%s_%s"%(j,systematic),MVA_Vars[systematic][j],tree))
1309     job.addpath('/MVAout')
1310     outfile = ROOT.TFile(job.getpath(), 'RECREATE')
1311     newtree = tree.CloneTree(0)
1312     #Setup Branches
1313     MVAbranches=[]
1314     for i in range(0,len(readers)):
1315     MVAbranches.append(array('f',[0]*9))
1316     newtree.Branch(MVAinfos[i].MVAname,MVAbranches[i],'nominal:JER_up:JER_down:JES_up:JES_down:beff_up:beff_down:bmis_up:bmis_down/F')
1317     print '\n--> ' + job.name +':'
1318     #progbar setup
1319     if nEntries >= longe:
1320     step=int(nEntries/longe)
1321     long=longe
1322     else:
1323     long=nEntries
1324     step = 1
1325     bar=progbar(long)
1326    
1327     #Fill event by event:
1328     for entry in range(0,nEntries):
1329     if entry % step == 0:
1330     bar.move()
1331     #load entry
1332     tree.GetEntry(entry)
1333     for systematic in systematics:
1334     for j in range(len( MVA_Vars['Nominal'])):
1335     MVA_var_buffer[j][0] = MVA_formulas[systematic][j].EvalInstance()
1336    
1337     for j in range(0,len(readers)):
1338     MVAbranches[j][systematics.index(systematic)] = readers[j].EvaluateMVA(MVAinfos[j].MVAname)
1339     #Fill:
1340     newtree.Fill()
1341    
1342     newtree.Write()
1343     newtree.Write()
1344     Count.Write()
1345     CountWithPU.Write()
1346     CountWithPU2011B.Write()
1347     outfile.Close()
1348    
1349     #DATA
1350     if job.type == 'DATA':
1351    
1352     #MVA Formulas
1353     MVA_formulas_Nominal = []
1354     #create TTreeFormulas
1355     for j in range(len( MVA_Vars['Nominal'])):
1356     MVA_formulas_Nominal.append(ROOT.TTreeFormula("MVA_formula%s_Nominal"%j, MVA_Vars['Nominal'][j],tree))
1357     job.addpath('/MVAout')
1358     outfile = ROOT.TFile(job.getpath(), 'RECREATE')
1359     newtree = tree.CloneTree(0)
1360     #Setup Branches
1361     MVAbranches=[]
1362     for i in range(0,len(readers)):
1363    
1364     MVAbranches.append(array('f',[0]))
1365     newtree.Branch(MVAinfos[i].MVAname,MVAbranches[i],'nominal/F')
1366     #progbar
1367     print '\n--> ' + job.name +':'
1368     if nEntries >= longe:
1369     step=int(nEntries/longe)
1370     long=longe
1371     else:
1372     long=nEntries
1373     step = 1
1374     bar=progbar(long)
1375    
1376     #Fill event by event:
1377     for entry in range(0,nEntries):
1378     if entry % step == 0:
1379     bar.move()
1380     #load entry
1381     tree.GetEntry(entry)
1382     #nominal:
1383     for j in range(len( MVA_Vars['Nominal'])):
1384     MVA_var_buffer[j][0] = MVA_formulas_Nominal[j].EvalInstance()
1385    
1386     for j in range(0,len(readers)):
1387     MVAbranches[j][0]= readers[j].EvaluateMVA(MVAinfos[j].MVAname)
1388    
1389    
1390     newtree.Fill()
1391     newtree.Write()
1392     newtree.Write()
1393     Count.Write()
1394     CountWithPU.Write()
1395     CountWithPU2011B.Write()
1396     outfile.Close()
1397    
1398     print '\n'
1399     infofile = open(Apath+'/MVAout/samples.info','w')
1400     pickle.dump(Ainfo,infofile)
1401     infofile.close()
1402    
1403    
1404     def Limit(path,var,data):
1405     print data
1406    
1407    
1408     plot=config.get('Limit',var)
1409    
1410     infofile = open(path+'/samples.info','r')
1411     info = pickle.load(infofile)
1412     infofile.close()
1413    
1414     options = plot.split(',')
1415    
1416    
1417     name=options[1]
1418     title = options[2]
1419     nBins=int(options[3])
1420     xMin=float(options[4])
1421     xMax=float(options[5])
1422    
1423     setup=config.get('Plot','setup')
1424     setup=setup.split(',')
1425    
1426    
1427     ROOToutname = options[6]
1428    
1429    
1430     outpath=config.get('Directories','limits')
1431    
1432    
1433     outfile = ROOT.TFile(outpath+ROOToutname+'.root', 'RECREATE')
1434     #Spuck out se Histograms for se Comination tool
1435     discr_names = ['Zudscg', 'Zbb', 'TTbar','VV', 'ST', 'Sig115', 'Wudscg', 'Wbb', 'QCD']
1436     data_name = ['data_obs']
1437    
1438    
1439    
1440    
1441     histos = []
1442     typs = []
1443     datas = []
1444     datatyps =[]
1445    
1446     for job in info:
1447     print job.name
1448     if job.type != 'DATA':
1449     print 'MC'
1450     hTemp, typ = getHistoFromTree2(job,options)
1451     histos.append(hTemp)
1452     typs.append(typ)
1453     elif job.name in data:
1454     print 'DATA'
1455     hTemp, typ = getHistoFromTree2(job,options)
1456     datas.append(hTemp)
1457     datatyps.append(typ)
1458    
1459    
1460    
1461    
1462     ROOT.gROOT.SetStyle("Plain")
1463     c = ROOT.TCanvas(name,title, 800, 600)
1464     MC_integral=0
1465     MC_entries=0
1466    
1467     for histo in histos:
1468     MC_integral+=histo.Integral()
1469     #MC_entries+=histo.GetEntries()
1470     print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
1471     #flow = MC_entries-MC_integral
1472     #if flow > 0:
1473     # print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
1474    
1475     #ORDER AND ADD TOGETHER
1476    
1477     ordnung=[]
1478     ordnungtyp=[]
1479     num=[0]*len(setup)
1480     for i in range(0,len(setup)):
1481     for j in range(0,len(histos)):
1482     if typs[j] == setup[i]:
1483     num[i]+=1
1484     ordnung.append(histos[j])
1485     ordnungtyp.append(typs[j])
1486    
1487     del histos
1488     del typs
1489    
1490     histos=ordnung
1491     typs=ordnungtyp
1492    
1493     for k in range(0,len(num)):
1494     for m in range(0,num[k]):
1495     if m > 0:
1496     histos[k].Add(histos[k+1],1)
1497     del histos[k+1]
1498     del typs[k+1]
1499    
1500    
1501    
1502     d1 = ROOT.TH1F('d1','d1',nBins,xMin,xMax)
1503    
1504     for i in range(0,len(datas)):
1505     d1.Add(datas[i],1)
1506     print "\033[1;32m\n\tDATA integral = %s\033[1;m"%d1.Integral()
1507     flow = d1.GetEntries()-d1.Integral()
1508     if flow > 0:
1509     print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
1510    
1511    
1512    
1513     for i in range(0,len(histos)):
1514     histos[i].SetName(discr_names[i])
1515     histos[i].SetDirectory(outfile)
1516     histos[i].Draw()
1517     print discr_names[i]
1518     print histos[i].Integral(0,nBins)
1519    
1520    
1521     #datas[0]: data_obs
1522     d1.SetName(data_name[0])
1523     d1.SetDirectory(outfile)
1524     print data_name[0]
1525     print d1.Integral(0,nBins)
1526     print d1.Integral()
1527     print d1.GetEntries()
1528    
1529     #write DATAcard
1530     f = open(outpath+'/vhbb_%s.txt'%ROOToutname,'w')
1531     f.write('imax\t1\tnumber of channels\njmax\t8\tnumber of backgrounds (\'*\' = automatic)\nkmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
1532     f.write('shapes * * %s.root $PROCESS $PROCESS$SYSTEMATIC\n\nbin\tZee\n\n'%ROOToutname)
1533     f.write('observation\t%s\n\n\n' %(d1.Integral()))
1534     f.write('bin\tZee\tZee\tZee\tZee\tZee\tZee\tZee\tZee\tZee\n')
1535     f.write('process\tSig115\tWudscg\tWbb\tZudscg\tZbb\tTTbar\tST\tVV\tQCD\n')
1536     f.write('process\t0\t1\t2\t3\t4\t5\t6\t7\t8\n')
1537     f.write('rate\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n\n' %(histos[5].Integral(),0,0,histos[0].Integral(),histos[1].Integral(),histos[2].Integral(),histos[4].Integral(),histos[3].Integral(),0)) #\t1.918\t0.000 0.000\t135.831 117.86 18.718 1.508\t7.015\t0.000
1538     f.write('lumi\tlnN\t1.045\t-\t-\t-\t-\t-\t1.045\t1.045\t1.045\npdf_qqbar\tlnN\t1.01\t-\t-\t-\t-\t-\t-\t1.01\t-\npdf_gg\tlnN\t-\t-\t-\t-\t-\t-\t1.01\t-\t1.01\nQCDscale_VH\tlnN\t1.04\t-\t-\t-\t-\t-\t-\t-\t-\nQCDscale_ttbar\tlnN\t-\t-\t-\t-\t-\t-\t1.06\t-\t-\nQCDscale_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.04\t-\nQCDscale_QCD\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t1.30\nCMS_vhbb_boost_EWK\tlnN\t1.05\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_boost_QCD\tlnN\t1.10\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_ST\tlnN\t-\t-\t-\t-\t-\t-\t1.29\t-\t-\nCMS_vhbb__VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.30\t-\nCMS_vhbb_WjLF_SF\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_WjHF_SF\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_ZjLF_SF\tlnN\t-\t-\t-\t1.06\t-\t-\t-\t-\t-\nCMS_vhbb_ZjHF_SF\tlnN\t-\t-\t-\t-\t1.17\t-\t-\t-\t-\nCMS_vhbb_TT_SF\tlnN\t-\t-\t-\t-\t-\t1.14\t-\t-\t-\nCMS_vhbb_QCD_SF\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_trigger_m\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_trigger_e\tlnN\t1.02\t-\t-\t-\t-\t-\t1.02\t1.02\t-\n')
1539     f.write('CMS_vhbb_trigger_MET\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_eff_m\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_eff_e\tlnN\t1.04\t-\t-\t-\t-\t-\t1.04\t1.04\t1.04\nCMS_toteff_b\tlnN\t1.10\t1.10\t1.00\t1.10\t1.00\t1.10\t1.10\t1.10\t1.10\nCMS_totscale_j\tlnN\t1.02\t-\t-\t-\t-\t-\t1.02\t1.02\t-\nCMS_totres_j\tlnN\t1.05\t1.03\t1.03\t1.03\t1.03\t1.03\t1.03\t1.05\t-\nCMS_vhbb_MET_nojets\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_VH_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_WjLF_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_WjHF_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhb_stats_ZjLF_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_ZjHF_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_TT_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_sT_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_VV_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_QCD_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1540     f.write('CMS_vhbb_stats_WjLF_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_WjHF_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhb_stats_ZjLF_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_ZjHF_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_TT_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_sT_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_VV_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_QCD_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_VH_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_WjLF_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_WjHF_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhb_stats_ZjLF_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_ZjHF_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_TT_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_sT_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_VV_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_QCD_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_VH_Zee\tlnN\t1.03\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_WjLF_Zee\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_WjHF_Zee\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhb_stats_ZjLF_Zee\tlnN\t-\t-\t-\t1.05\t-\t-\t-\t-\t-\nCMS_vhbb_stats_ZjHF_Zee\tlnN\t-\t-\t-\t-\t1.07\t-\t-\t-\t-\nCMS_vhbb_stats_TT_Zee\tlnN\t-\t-\t-\t-\t-\t1.06\t-\t-\t-\nCMS_vhbb_stats_sT_Zee\tlnN\t-\t-\t-\t-\t-\t-\t1.30\t-\t-\nCMS_vhbb_stats_Diboson_Zee\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.06\t-\nCMS_vhbb_stats_QCD_Zee\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_VH_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_WjLF_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_WjHF_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhb_stats_ZjLF_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_ZjHF_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_TT_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_sT_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_VV_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\nCMS_vhbb_stats_QCD_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1541     f.close()
1542    
1543     #dunnmies
1544     dummies=[]
1545     #Wlight,Wbb,QCD
1546     for i in range(6,9):
1547     dummy = ROOT.TH1F(discr_names[i], "discriminator", nBins, xMin, xMax)
1548     dummy.SetDirectory(outfile)
1549     dummy.Draw()
1550     dummies.append(dummy)
1551     print discr_names[i]
1552     #Wbb
1553     #dummy = ROOT.TH1F(discr_names[7], "discriminator", div, discrMin, discrMax)
1554     #dummies.append(dummy)
1555     #QCD
1556     #dummy = ROOT.TH1F(discr_names[8], "discriminator", div, discrMin, discrMax)
1557     #dummies.append(dummy)
1558    
1559     #Write to file
1560     outfile.Write()
1561     outfile.Close()
1562    
1563    
1564     def writeWorkspace(path,var,data):
1565    
1566    
1567     plot=config.get('Limit',var)
1568    
1569     infofile = open(path+'/samples.info','r')
1570     info = pickle.load(infofile)
1571     infofile.close()
1572    
1573     options = plot.split(',')
1574    
1575    
1576     name=options[1]
1577     title = options[2]
1578     nBins=int(options[3])
1579     xMin=float(options[4])
1580     xMax=float(options[5])
1581    
1582     setup=config.get('Plot','setup')
1583     setup=setup.split(',')
1584    
1585    
1586     ROOToutname = options[6]
1587    
1588    
1589     outpath=config.get('Directories','limits')
1590    
1591    
1592     outfile = ROOT.TFile(outpath+ROOToutname+'_WS.root', 'RECREATE')
1593     discr_names = ['Zudscg', 'Zbb', 'TTbar','VV', 'ST', 'Sig115', 'Wudscg', 'Wbb', 'QCD']
1594     data_name = ['data_obs']
1595    
1596     WS = ROOT.RooWorkspace('Zee','Zee')
1597     print 'WS initialized'
1598    
1599     disc= ROOT.RooRealVar('BDT','BDT',-1,1)
1600     obs = ROOT.RooArgList(disc)
1601    
1602    
1603     histos = []
1604     typs = []
1605     datas = []
1606     datatyps =[]
1607    
1608     for job in info:
1609     #print job.name
1610     if job.type != 'DATA':
1611     #print 'MC'
1612     hTemp, typ = getHistoFromTree2(job,options)
1613     histos.append(hTemp)
1614     typs.append(typ)
1615     elif job.name in data:
1616     #print 'DATA'
1617     hTemp, typ = getHistoFromTree2(job,options)
1618     datas.append(hTemp)
1619     datatyps.append(typ)
1620    
1621    
1622    
1623    
1624     ROOT.gROOT.SetStyle("Plain")
1625     c = ROOT.TCanvas(name,title, 800, 600)
1626     MC_integral=0
1627     MC_entries=0
1628    
1629     for histo in histos:
1630     MC_integral+=histo.Integral()
1631     #MC_entries+=histo.GetEntries()
1632     print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
1633     #flow = MC_entries-MC_integral
1634     #if flow > 0:
1635     # print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
1636    
1637     #ORDER AND ADD TOGETHER
1638    
1639     ordnung=[]
1640     ordnungtyp=[]
1641     num=[0]*len(setup)
1642     for i in range(0,len(setup)):
1643     for j in range(0,len(histos)):
1644     if typs[j] == setup[i]:
1645     num[i]+=1
1646     ordnung.append(histos[j])
1647     ordnungtyp.append(typs[j])
1648    
1649     del histos
1650     del typs
1651    
1652     histos=ordnung
1653     typs=ordnungtyp
1654    
1655     for k in range(0,len(num)):
1656     for m in range(0,num[k]):
1657     if m > 0:
1658     histos[k].Add(histos[k+1],1)
1659     del histos[k+1]
1660     del typs[k+1]
1661    
1662    
1663    
1664     d1 = ROOT.TH1F('d1','d1',nBins,xMin,xMax)
1665    
1666     for i in range(0,len(datas)):
1667     d1.Add(datas[i],1)
1668     print "\033[1;32m\n\tDATA integral = %s\033[1;m"%d1.Integral()
1669     flow = d1.GetEntries()-d1.Integral()
1670     if flow > 0:
1671     print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
1672    
1673    
1674    
1675    
1676    
1677    
1678     for i in range(0,len(histos)):
1679     histos[i].SetName(discr_names[i])
1680     histos[i].SetDirectory(outfile)
1681     histos[i].Draw()
1682    
1683    
1684    
1685    
1686    
1687     statUp = histos[i].Clone()
1688     statDown = histos[i].Clone()
1689     #shift up and down with statistical error
1690     for j in range(histos[i].GetNbinsX()):
1691     statUp.SetBinContent(j,statUp.GetBinContent(j)+statUp.GetBinError(j))
1692     statDown.SetBinContent(j,statDown.GetBinContent(j)-statDown.GetBinError(j))
1693     statUp.SetName('%sStatsUp'%discr_names[i])
1694     statDown.SetName('%sStatsDown'%discr_names[i])
1695    
1696    
1697     histPdf = ROOT.RooDataHist(discr_names[i],discr_names[i],obs,histos[i])
1698    
1699     #UP stats of MCs
1700     RooStatsUp = ROOT.RooDataHist('%sStatsUp'%discr_names[i],'%sStatsUp'%discr_names[i],obs, statUp)
1701     #DOWN stats of MCs
1702     RooStatsDown = ROOT.RooDataHist('%sStatsDown'%discr_names[i],'%sStatsDown'%discr_names[i],obs, statDown)
1703    
1704    
1705     getattr(WS,'import')(histPdf)
1706     getattr(WS,'import')(RooStatsUp)
1707     getattr(WS,'import')(RooStatsDown)
1708    
1709    
1710    
1711     frame=disc.frame()
1712    
1713    
1714     ROOT.RooAbsData.plotOn(histPdf,frame)
1715     frame.Draw()
1716    
1717     c.Print('~/Hbb/WStest/%s.png'%discr_names[i])
1718    
1719    
1720    
1721    
1722     #print discr_names[i]
1723     #print histos[i].Integral(0,nBins)
1724    
1725    
1726     #datas[0]: data_obs
1727     d1.SetName(data_name[0])
1728     d1.SetDirectory(outfile)
1729     #print data_name[0]
1730     #print d1.Integral(0,nBins)
1731     #print d1.Integral()
1732     #print d1.GetEntries()
1733    
1734     #write DATAcard
1735     f = open(outpath+'/vhbb_%s_WS.txt'%ROOToutname,'w')
1736     f.write('imax\t1\tnumber of channels\n')
1737     f.write('jmax\t8\tnumber of backgrounds (\'*\' = automatic)\n')
1738     f.write('kmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
1739    
1740     f.write('shapes * * %s_WS.root $CHANNEL:$PROCESS $CHANNEL:$PROCESS$SYSTEMATIC\n\n'%ROOToutname)
1741     f.write('bin\tZee\n\n')
1742     f.write('observation\t%s\n\n'%d1.Integral())
1743     f.write('bin\tZee\tZee\tZee\tZee\tZee\tZee\tZee\tZee\tZee\n')
1744     f.write('process\tSig115\tWudscg\tWbb\tZudscg\tZbb\tTTbar\tST\tVV\tQCD\n')
1745     f.write('process\t0\t1\t2\t3\t4\t5\t6\t7\t8\n')
1746     f.write('rate\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n'%(histos[5].Integral(),0,0,histos[0].Integral(),histos[1].Integral(),histos[2].Integral(),histos[4].Integral(),histos[3].Integral(),0)) #\t1.918\t0.000 0.000\t135.831 117.86 18.718 1.508\t7.015\t0.000
1747     f.write('lumi\tlnN\t1.045\t-\t-\t-\t-\t-\t1.045\t1.045\t1.045\n\n')
1748     f.write('pdf_qqbar\tlnN\t1.01\t-\t-\t-\t-\t-\t-\t1.01\t-\n')
1749     f.write('pdf_gg\tlnN\t-\t-\t-\t-\t-\t-\t1.01\t-\t1.01\n')
1750     f.write('QCDscale_VH\tlnN\t1.04\t-\t-\t-\t-\t-\t-\t-\t-\n')
1751     f.write('QCDscale_ttbar\tlnN\t-\t-\t-\t-\t-\t-\t1.06\t-\t-\n')
1752     f.write('QCDscale_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.04\t-\n')
1753     f.write('QCDscale_QCD\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t1.30\n')
1754     f.write('CMS_vhbb_boost_EWK\tlnN\t1.05\t-\t-\t-\t-\t-\t-\t-\t-\n')
1755     f.write('CMS_vhbb_boost_QCD\tlnN\t1.10\t-\t-\t-\t-\t-\t-\t-\t-\n')
1756     f.write('CMS_vhbb_ST\tlnN\t-\t-\t-\t-\t-\t-\t1.29\t-\t-\n')
1757     f.write('CMS_vhbb__VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.30\t-\n')
1758     f.write('CMS_vhbb_WjLF_SF\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1759     f.write('CMS_vhbb_WjHF_SF\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1760     f.write('CMS_vhbb_ZjLF_SF\tlnN\t-\t-\t-\t1.06\t-\t-\t-\t-\t-\n')
1761     f.write('CMS_vhbb_ZjHF_SF\tlnN\t-\t-\t-\t-\t1.17\t-\t-\t-\t-\n')
1762     f.write('CMS_vhbb_TT_SF\tlnN\t-\t-\t-\t-\t-\t1.14\t-\t-\t-\n')
1763     f.write('CMS_vhbb_QCD_SF\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1764     f.write('CMS_trigger_m\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1765     f.write('CMS_trigger_e\tlnN\t1.02\t-\t-\t-\t-\t-\t1.02\t1.02\t-\n')
1766     f.write('CMS_vhbb_trigger_MET\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1767     f.write('CMS_eff_m\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1768     f.write('CMS_eff_e\tlnN\t1.04\t-\t-\t-\t-\t-\t1.04\t1.04\t1.04\n')
1769     f.write('CMS_toteff_b\tlnN\t1.10\t1.10\t1.00\t1.10\t1.00\t1.10\t1.10\t1.10\t1.10\n')
1770     f.write('CMS_totscale_j\tlnN\t1.02\t-\t-\t-\t-\t-\t1.02\t1.02\t-\n')
1771     f.write('CMS_totres_j\tlnN\t1.05\t1.03\t1.03\t1.03\t1.03\t1.03\t1.03\t1.05\t-\n')
1772     f.write('CMS_vhbb_MET_nojets\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1773     f.write('CMS_vhbb_stats_VH_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1774     f.write('CMS_vhbb_stats_WjLF_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1775     f.write('CMS_vhbb_stats_WjHF_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1776     f.write('CMS_vhb_stats_ZjLF_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1777     f.write('CMS_vhbb_stats_ZjHF_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1778     f.write('CMS_vhbb_stats_TT_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1779     f.write('CMS_vhbb_stats_sT_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1780     f.write('CMS_vhbb_stats_VV_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1781     f.write('CMS_vhbb_stats_QCD_Wmn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1782     f.write('CMS_vhbb_stats_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1783     f.write('CMS_vhbb_stats_WjLF_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1784     f.write('CMS_vhbb_stats_WjHF_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1785     f.write('CMS_vhb_stats_ZjLF_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1786     f.write('CMS_vhbb_stats_ZjHF_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1787     f.write('CMS_vhbb_stats_TT_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1788     f.write('CMS_vhbb_stats_sT_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1789     f.write('CMS_vhbb_stats_VV_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1790     f.write('CMS_vhbb_stats_QCD_Wen\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1791     f.write('CMS_vhbb_stats_VH_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1792     f.write('CMS_vhbb_stats_WjLF_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1793     f.write('CMS_vhbb_stats_WjHF_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1794     f.write('CMS_vhb_stats_ZjLF_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1795     f.write('CMS_vhbb_stats_ZjHF_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1796     f.write('CMS_vhbb_stats_TT_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1797     f.write('CMS_vhbb_stats_sT_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1798     f.write('CMS_vhbb_stats_VV_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1799     f.write('CMS_vhbb_stats_QCD_Zmm\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1800     f.write('CMS_vhbb_stats_VH_Zee\tlnN\t1.03\t-\t-\t-\t-\t-\t-\t-\t-\n')
1801     f.write('CMS_vhbb_stats_WjLF_Zee\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1802     f.write('CMS_vhbb_stats_WjHF_Zee\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1803     f.write('CMS_vhb_stats_ZjLF_Zee\tlnN\t-\t-\t-\t1.05\t-\t-\t-\t-\t-\n')
1804     f.write('CMS_vhbb_stats_ZjHF_Zee\tlnN\t-\t-\t-\t-\t1.07\t-\t-\t-\t-\n')
1805     f.write('CMS_vhbb_stats_TT_Zee\tlnN\t-\t-\t-\t-\t-\t1.06\t-\t-\t-\n')
1806     f.write('CMS_vhbb_stats_sT_Zee\tlnN\t-\t-\t-\t-\t-\t-\t1.30\t-\t-\n')
1807     f.write('CMS_vhbb_stats_Diboson_Zee\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.06\t-\n')
1808     f.write('CMS_vhbb_stats_QCD_Zee\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1809     f.write('CMS_vhbb_stats_VH_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1810     f.write('CMS_vhbb_stats_WjLF_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1811     f.write('CMS_vhbb_stats_WjHF_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1812     f.write('CMS_vhb_stats_ZjLF_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1813     f.write('CMS_vhbb_stats_ZjHF_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1814     f.write('CMS_vhbb_stats_TT_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1815     f.write('CMS_vhbb_stats_sT_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1816     f.write('CMS_vhbb_stats_VV_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1817     f.write('CMS_vhbb_stats_QCD_Znn\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
1818    
1819     f.write('Stats\tshape\t1.0\t-\t-\t-\t-\t-\t-\t-\t-\n')
1820     f.write('Stats\tshape\t-\t1.0\t-\t-\t-\t-\t-\t-\t-\n')
1821     f.write('Stats\tshape\t-\t-\t1.0\t-\t-\t-\t-\t-\t-\n')
1822     f.write('Stats\tshape\t-\t-\t-\t1.0\t-\t-\t-\t-\t-\n')
1823     f.write('Stats\tshape\t-\t-\t-\t-\t1.0\t-\t-\t-\t-\n')
1824     f.write('Stats\tshape\t-\t-\t-\t-\t-\t1.0\t-\t-\t-\n')
1825     f.write('Stats\tshape\t-\t-\t-\t-\t-\t-\t1.0\t-\t-\n')
1826     f.write('Stats\tshape\t-\t-\t-\t-\t-\t-\t-\t1.0\t-\n')
1827     f.write('Stats\tshape\t-\t-\t-\t-\t-\t-\t-\t-\t1.0\n')
1828    
1829    
1830     f.close()
1831    
1832     #dunnmies
1833     #Wlight,Wbb,QCD
1834     for i in range(6,9):
1835     dummy = ROOT.TH1F(discr_names[i], "discriminator", nBins, xMin, xMax)
1836     dummy.SetDirectory(outfile)
1837     dummy.Draw()
1838    
1839     #nominal
1840     histPdf = ROOT.RooDataHist(discr_names[i],discr_names[i],obs,dummy)
1841     #UP stats of MCs
1842     RooStatsUp = ROOT.RooDataHist('%sStatsUp'%discr_names[i],'%sStatsUp'%discr_names[i],obs, dummy)
1843     #DOWN stats of MCs
1844     RooStatsDown = ROOT.RooDataHist('%sStatsDown'%discr_names[i],'%sStatsDown'%discr_names[i],obs, dummy)
1845    
1846    
1847     getattr(WS,'import')(histPdf)
1848     getattr(WS,'import')(RooStatsUp)
1849     getattr(WS,'import')(RooStatsDown)
1850     #print discr_names[i]
1851     #Wbb
1852     #dummy = ROOT.TH1F(discr_names[7], "discriminator", div, discrMin, discrMax)
1853     #dummies.append(dummy)
1854     #QCD
1855     #dummy = ROOT.TH1F(discr_names[8], "discriminator", div, discrMin, discrMax)
1856     #dummies.append(dummy)
1857    
1858     #Write to file
1859    
1860    
1861    
1862    
1863     #HISTOGRAMM of DATA
1864     #ROOT.RooDataHist('data_obsHist','',RooArgList,??)
1865     histPdf = ROOT.RooDataHist('data_obs','data_obs',obs,d1)
1866     ROOT.RooAbsData.plotOn(histPdf,frame)
1867     frame.Draw()
1868    
1869     c.Print('~/Hbb/WStest/d1.png')
1870     #IMPORT
1871     getattr(WS,'import')(histPdf)
1872    
1873     #Number of Obs?
1874     nObs = int(d1.Integral())
1875    
1876    
1877    
1878    
1879     '''
1880    
1881     theStatsUp = []
1882     theStatsDown = []
1883    
1884     #LOOP over MCsamples BKG
1885     for i in range(0,len(self.__theStacks)):
1886     #name = 'ZjLF'
1887     name = '%s%s' %(self.__dcRepMap[self.__datasets[i]],self.__writeCombination) #what is self writeCombination??, assume ''
1888     print name
1889     self.__theStacks[i].SetName(name)
1890    
1891     #HISTOGRAMM of MCs
1892     #ROOT.RooDataHist('ZjLF','',RooArgList,??)
1893     histPdf = ROOT.RooDataHist(name,'',obs,self.__theStacks[i])
1894     #IMPORT
1895     getattr(self.__w,'import')(histPdf)
1896    
1897    
1898     if self.__writeCombination == '':
1899     self.__dcRepMap['n%s'%(self.__datasets[i])] = self.__theStacks[i].Integral()
1900     name = '%s%s%s%s' %(self.__dcRepMap[self.__datasets[i]],'_CMS_vhbb_stats_',self.__dcRepMap[self.__datasets[i]],'_%(bin)s'%self.__dcRepMap)
1901     statUp = self.__theStacks[i].Clone()
1902     statDown = self.__theStacks[i].Clone()
1903     #shift up and down with statistical error
1904     for j in range(self.__theStacks[i].GetNbinsX()):
1905     statUp.SetBinContent(j,statUp.GetBinContent(j)+statUp.GetBinError(j))
1906     statDown.SetBinContent(j,statDown.GetBinContent(j)-statDown.GetBinError(j))
1907     theStatsUp.append(statUp)
1908     theStatsDown.append(statDown)
1909     theStatsUp[i].SetName('%s%s' %(name,'Up'))
1910     theStatsDown[i].SetName('%s%s' %(name,'Down'))
1911     #UP stats of MCs
1912     theRooStatsUp = ROOT.RooDataHist('%s%s' %(name,'Up'),'',obs, theStatsUp[i])
1913     #DOWN stats of MCs
1914     theRooStatsDown = ROOT.RooDataHist('%s%s' %(name,'Down'),'',obs, theStatsDown[i])
1915     getattr(self.__w,'import')(theRooStatsUp)
1916     getattr(self.__w,'import')(theRooStatsDown)
1917    
1918    
1919     #overlays=signal SIG
1920     #OVERLAYS??
1921     theOStatsUp = []
1922     theOStatsDown = []
1923     for i in range(0,len(self.__theOverlays)):
1924     name = '%s%s' %(self.__dcRepMap[self.__overlays[i]],self.__writeCombination)
1925     self.__theOverlays[i].SetName(name)
1926     histPdf = ROOT.RooDataHist(name,'',obs,self.__theOverlays[i])
1927     getattr(self.__w,'import')(histPdf)
1928     if self.__writeCombination == '':
1929     self.__dcRepMap['nSig'] = self.__theOverlays[i].Integral()
1930     #e.g. name=TTbar_CMS_vhbb_stats_TTbar_ZeeUp
1931     name = '%s%s%s%s' %(self.__dcRepMap[self.__overlays[i]],'_CMS_vhbb_stats_',self.__dcRepMap[self.__overlays[i]],'_%(bin)s'%self.__dcRepMap)
1932     statUp = self.__theOverlays[i].Clone()
1933     statDown = self.__theOverlays[i].Clone()
1934     for j in range(self.__theOverlays[i].GetNbinsX()):
1935     statUp.SetBinContent(j,statUp.GetBinContent(j)+statUp.GetBinError(j))
1936     statDown.SetBinContent(j,statDown.GetBinContent(j)-statDown.GetBinError(j))
1937     theOStatsUp.append(statUp)
1938     theOStatsDown.append(statDown)
1939     theOStatsUp[i].SetName('%s%s' %(name,'Up'))
1940     theOStatsDown[i].SetName('%s%s' %(name,'Down'))
1941     theRooStatsUp = ROOT.RooDataHist('%s%s' %(name,'Up'),'',obs, theOStatsUp[i])
1942     theRooStatsDown = ROOT.RooDataHist('%s%s' %(name,'Down'),'',obs, theOStatsDown[i])
1943     getattr(self.__w,'import')(theRooStatsUp)
1944     getattr(self.__w,'import')(theRooStatsDown)
1945    
1946     '''
1947    
1948     WS.writeToFile(outpath+ROOToutname+'_WS.root')
1949     #WS.writeToFile("testWS.root")
1950    
1951    
1952    
1953    
1954    
1955    
1956    
1957    
1958     def SysPlot(mode,systematic):
1959    
1960     ROOT.gROOT.SetStyle("Plain")
1961     c = ROOT.TCanvas('title','title', 800, 600)
1962     ROOT.gPad.SetTicks(1,1)
1963    
1964    
1965     #systematic='JER'
1966    
1967     #if mode == 'test':
1968     # type = 'TMVAClassification_nov10BDTCatnaJet3_shuffled'
1969    
1970     #if mode == 'test2':
1971     # type = 'TMVAClassification_nov10BDT_shuffled'
1972    
1973     print 'ok, i plot the MVA output for you...'
1974     #namehisto = 'taskTMVAClassification_BDTCatnaJet3loose'
1975     namehisto = task+type
1976     rebin = 100
1977     if mode == 'test': path=treePath+'/test'
1978     if mode == 'Top': path=treePath+'/Top'
1979     if mode == 'Zlight': path=treePath+'/Zlight'
1980     if mode == 'Zbb': path=treePath+'/Zbb'
1981     if mode == 'Signal': path=treePath+'/Signal'
1982    
1983     MVAtitle=mode
1984     nBins=div/rebin
1985    
1986     Ntotal = ROOT.TH1F(systematic,systematic,nBins,discrMin,discrMax)
1987     Utotal = ROOT.TH1F('Utotal','Utotal',nBins,discrMin,discrMax)
1988     Dtotal = ROOT.TH1F('Dtotal','Dtotal',nBins,discrMin,discrMax)
1989    
1990     for job in jobs: #jobs:
1991     jobN= path +'/MVA_'+training+'_'+MVAtitle+'.' + job +'.root'
1992     jobU= path +'/MVA_'+training+'_'+MVAtitle+'.' + job +'.'+systematic+'_up.root'
1993     jobD= path +'/MVA_'+training+'_'+MVAtitle+'.' + job +'.'+systematic+'_down.root'
1994     print jobN
1995     l = ROOT.TLegend(0.28, 0.73, 0.38, 0.88)
1996     #hTemp = getHistoFromTree(path,job2,0)
1997    
1998     N = ROOT.TFile(jobN, 'OPEN')
1999     NHist = N.Get(namehisto)
2000     NHist.Rebin(rebin)
2001     NHist.SetDirectory(0)
2002     NHist.SetLineColor(1)
2003     NHist.SetMarkerStyle(8)
2004     NHist.SetStats(0)
2005     NHist.SetTitle('MVA '+systematic+' '+ legenden[jobs.index(job)])
2006     Ntotal.Add(NHist)
2007     l.AddEntry(NHist,'nominal','PL')
2008    
2009     U = ROOT.TFile(jobU, 'OPEN')
2010     UHist = U.Get(namehisto)
2011     UHist.Rebin(rebin)
2012     UHist.SetDirectory(0)
2013     UHist.SetLineColor(4)
2014     UHist.SetLineStyle(4)
2015     UHist.SetLineWidth(2)
2016     l.AddEntry(UHist,'up','PL')
2017     Utotal.Add(UHist)
2018    
2019     D = ROOT.TFile(jobD, 'OPEN')
2020     DHist = D.Get(namehisto)
2021     DHist.Rebin(rebin)
2022     DHist.SetDirectory(0)
2023     DHist.SetLineColor(2)
2024     DHist.SetLineStyle(3)
2025     DHist.SetLineWidth(2)
2026     l.AddEntry(DHist,'down','PL')
2027     Dtotal.Add(DHist)
2028    
2029     NHist.Draw("P0")
2030     NHist.Draw("same")
2031     UHist.Draw("same")
2032     DHist.Draw("same")
2033     l.SetFillColor(0)
2034     l.SetBorderSize(0)
2035     l.Draw()
2036     title= mode + type + legenden[jobs.index(job)] +systematic
2037     name = '%s/Stack/%s.png' %(plotPath,title)
2038     c.Print(name)
2039     N.Close()
2040     U.Close()
2041     D.Close()
2042    
2043     Ntotal.SetMarkerStyle(8)
2044     Ntotal.SetLineColor(1)
2045     Ntotal.SetStats(0)
2046     Ntotal.Draw("P0")
2047     Ntotal.Draw("same")
2048     Utotal.SetLineColor(4)
2049     Utotal.SetLineStyle(4)
2050     Utotal.SetLineWidth(2)
2051     Utotal.Draw("same")
2052     Dtotal.SetLineColor(2)
2053     Dtotal.SetLineStyle(3)
2054     Dtotal.Draw("same")
2055     Dtotal.SetLineWidth(2)
2056     l.Draw()
2057    
2058     title= mode + type +systematic
2059     name = '%s/Stack/%s.png' %(plotPath,title)
2060     c.Print(name)
2061    
2062     def newFoM(path,var):
2063    
2064    
2065    
2066    
2067    
2068     plot=config.get('FoM',var)
2069    
2070     infofile = open(path+'/samples.info','r')
2071     info = pickle.load(infofile)
2072     infofile.close()
2073    
2074    
2075    
2076     options = plot.split(',')
2077     name=options[1]
2078     title = options[2]
2079     nBins=int(options[3])
2080     xMin=float(options[4])
2081     xMax=float(options[5])
2082    
2083     bkgs=config.get('FoM','BKG')
2084     bkgs=bkgs.split(' ')
2085    
2086     sigs=config.get('FoM','SIG')
2087     sigs=sigs.split(' ')
2088    
2089    
2090     ROOT.gROOT.SetStyle("Plain")
2091     c = ROOT.TCanvas(title,title, 800, 600)
2092     ROOT.gPad.SetTicks(1,1)
2093    
2094    
2095    
2096    
2097     print '\nProducing Plot of %s\n'%title
2098    
2099    
2100     histos = []
2101    
2102     for job in info:
2103     if job.name in bkgs:
2104     hTemp, typ = getHistoFromTree2(job,options)
2105     histos.append(hTemp)
2106    
2107     for job in info:
2108     if job.name in sigs:
2109     hTemp, typ = getHistoFromTree2(job,options)
2110     histos.append(hTemp)
2111    
2112     for i in range(1,len(bkgs)):
2113     histos[0].Add(histos[1],1)
2114     del histos[1]
2115    
2116    
2117     for i in range(1,len(sigs)):
2118     histos[1].Add(histos[2],1)
2119     del histos[2]
2120    
2121    
2122    
2123     fig = []
2124     C=[]
2125     print '\n\t--> Info:'
2126     B=histos[0].Integral()
2127     print '\t Background count is %s' %B
2128     S=histos[1].Integral()
2129     print '\t Signal count is %s\n' %S
2130     F=[]
2131     print 'nbins %s' %nBins
2132     print 'size %s' %histos[0].GetSize()
2133     for i in range(0,nBins):
2134     #print S
2135     #print B
2136     if B >= 0:
2137     FOM = S/(1.5+sqrt(B)+0.2*B)
2138     F.append(FOM)
2139     #print (S/(1.5+sqrt(B)+0.2*B))
2140     print 'S = %s, B = %s, FoM = %s' %(S,B,FOM)
2141     C.append(histos[0].GetBinCenter(i+1))
2142     else: print 'S %s, B %s' %(S,B)
2143     B=B-histos[0].GetBinContent(i+1)
2144     S=S-histos[1].GetBinContent(i+1)
2145    
2146     x = array('f', C)
2147     y = array('f', F)
2148     gr1 = ROOT.TGraph(len(x), x,y)
2149     gr1.SetTitle(title)
2150     gr1.Draw('APL')
2151     gr1.GetXaxis().SetTitle('BDT Cut')
2152     gr1.GetYaxis().SetTitle('FoM')
2153     gr1.GetXaxis().SetRangeUser(xMin,xMax)
2154     name = '%s/%s' %(config.get('Directories','plotpath'),options[6])
2155     c.Print(name)
2156    
2157     #*********************Actually DO STH*******************************
2158     #getList(signalFiles)
2159     if sys.argv[1] == 'limitWS': writeWorkspace(sys.argv[2],sys.argv[3],sys.argv[4])
2160     if sys.argv[1] == 'newtrain': newTraining(sys.argv[2],sys.argv[3])
2161     if sys.argv[1] == 'eval': evaluate(sys.argv[2],sys.argv[3])#,sys.argv[4])
2162     if sys.argv[1] == 'evalMulti': evalMulti(sys.argv[2],sys.argv[3])
2163     if sys.argv[1] == 'limit': Limit(sys.argv[2],sys.argv[3],sys.argv[4])
2164     if sys.argv[1] == 'plot': plot()
2165     if sys.argv[1] == 'allplots': allplots()
2166     if sys.argv[1] == 'comp': createComparison()
2167     if sys.argv[1] == 'compare': treeCompare(sys.argv[2],sys.argv[3])
2168    
2169     if sys.argv[1] == 'MVAplot': MVAstack(sys.argv[2])
2170     if sys.argv[1] == 'SysPlot': SysPlot(sys.argv[2],sys.argv[3])
2171     if sys.argv[1] == 'copy': CutCopy(sys.argv[2])
2172     if sys.argv[1] == 'FoM': newFoM(sys.argv[2],sys.argv[3])
2173     if sys.argv[1] == 'shuffle': shuffle()
2174     if sys.argv[1] == 'SuperShuffle': SuperShuffle()
2175     if sys.argv[1] == 'addsys': AddSystematics(sys.argv[2])
2176     if sys.argv[1] == 'addcut': Addcut(sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5])
2177     if sys.argv[1] == 'addsinglecut': Addsinglecut(sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5])
2178     if sys.argv[1] == 'addfile': AddFile(sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5],sys.argv[6])
2179     if sys.argv[1] == 'stack': treeStack(sys.argv[2],sys.argv[3],sys.argv[4])