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

# Content
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])