ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/workspace_datacard.py
Revision: 1.47
Committed: Sun Oct 14 09:53:49 2012 UTC (12 years, 7 months ago) by peller
Content type: text/x-python
Branch: MAIN
Changes since 1.46: +6 -1 lines
Log Message:
Mjj

File Contents

# User Rev Content
1 peller 1.8 #!/usr/bin/env python
2 peller 1.1 import os
3 peller 1.11 import sys
4 peller 1.1 import ROOT
5     from ROOT import TFile
6     from array import array
7     from math import sqrt
8     from copy import copy
9     #suppres the EvalInstace conversion warning bug
10     import warnings
11     warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
12 nmohr 1.31 from optparse import OptionParser
13 nmohr 1.20 from BetterConfigParser import BetterConfigParser
14 peller 1.1 from samplesclass import sample
15     from mvainfos import mvainfo
16     import pickle
17     from progbar import progbar
18     from printcolor import printc
19 peller 1.13 from gethistofromtree import getHistoFromTree, orderandadd
20 peller 1.1
21     #CONFIGURE
22 nmohr 1.31 argv = sys.argv
23     parser = OptionParser()
24     parser.add_option("-P", "--path", dest="path", default="",
25     help="path to samples")
26     parser.add_option("-V", "--var", dest="variable", default="",
27     help="variable for shape analysis")
28     parser.add_option("-C", "--config", dest="config", default=[], action="append",
29     help="configuration file")
30     (opts, args) = parser.parse_args(argv)
31     if opts.config =="":
32     opts.config = "config"
33     print opts.config
34 nmohr 1.20 config = BetterConfigParser()
35 nmohr 1.31 config.read(opts.config)
36     anaTag = config.get("Analysis","tag")
37 peller 1.30
38 peller 1.36
39     # -------------------- parsing configuration and options: (an ugly spaghetti code section) ----------------------------------------------------------------------
40 peller 1.1 #get locations:
41     Wdir=config.get('Directories','Wdir')
42 peller 1.42 vhbbpath=config.get('Directories','vhbbpath')
43 nmohr 1.46 samplesinfo=config.get('Directories','samplesinfo')
44 peller 1.1 #systematics
45     systematics=config.get('systematics','systematics')
46     systematics=systematics.split(' ')
47     weightF=config.get('Weights','weightF')
48 nmohr 1.31 path=opts.path
49     var=opts.variable
50 peller 1.1 plot=config.get('Limit',var)
51 nmohr 1.46 infofile = open(samplesinfo,'r')
52 peller 1.1 info = pickle.load(infofile)
53     infofile.close()
54     options = plot.split(',')
55 bortigno 1.7 if len(options) < 12:
56     print "You have to choose option[11]: either Mjj or BDT"
57     sys.exit("You have to choose option[11]: either Mjj or BDT")
58 peller 1.1 name=options[1]
59     title = options[2]
60     nBins=int(options[3])
61     xMin=float(options[4])
62     xMax=float(options[5])
63 peller 1.29 SIG=options[9]
64 peller 1.1 data=options[10]
65 bortigno 1.7 anType=options[11]
66 peller 1.13 RCut=options[7]
67 peller 1.29 setup=eval(config.get('LimitGeneral','setup'))
68 peller 1.13 ROOToutname = options[6]
69     outpath=config.get('Directories','limits')
70     outfile = ROOT.TFile(outpath+'vhbb_TH_'+ROOToutname+'.root', 'RECREATE')
71 peller 1.29 systematicsnaming=eval(config.get('LimitGeneral','systematicsnaming7TeV'))
72 peller 1.44
73     TrainFlag = eval(config.get('Analysis','TrainFlag'))
74    
75    
76 peller 1.30 if anaTag =='8TeV':
77 peller 1.29 systematicsnaming=eval(config.get('LimitGeneral','systematicsnaming8TeV'))
78 peller 1.44 elif not anaTag =='7TeV':
79 peller 1.30 print "What is your Analysis Tag in config? (anaTag)"
80     sys.exit("What is your Analysis Tag in config? (anaTag)")
81 peller 1.29 scaling=eval(config.get('LimitGeneral','scaling'))
82 peller 1.44
83     if TrainFlag:
84     MC_rescale_factor=2.
85     print 'I RESCALE BY 2.0'
86     else: MC_rescale_factor = 1.
87    
88 peller 1.29 rescaleSqrtN=eval(config.get('LimitGeneral','rescaleSqrtN'))
89 peller 1.11 if 'RTight' in RCut:
90 nmohr 1.17 Datacradbin=options[10]
91 peller 1.11 elif 'RMed' in RCut:
92 nmohr 1.17 Datacradbin=options[10]
93 peller 1.11 else:
94     Datacradbin=options[10]
95 peller 1.36 blind=eval(config.get('LimitGeneral','blind'))
96     BKGlist = eval(config.get('LimitGeneral','BKG'))
97     #Groups for adding samples together
98     Group = eval(config.get('LimitGeneral','Group'))
99     #naming for DC
100     Dict= eval(config.get('LimitGeneral','Dict'))
101     weightF_sys = eval(config.get('LimitGeneral','weightF_sys'))
102     binstat = eval(config.get('LimitGeneral','binstat'))
103 nmohr 1.37 addSample_sys = None if not config.has_option('LimitGeneral','addSample_sys') else eval(config.get('LimitGeneral','addSample_sys'))
104 peller 1.36 bdt = False
105     mjj = False
106     #print str(anType)
107     #print len(options)
108     if str(anType) == 'BDT':
109     bdt = True
110     systematics = eval(config.get('LimitGeneral','sys_BDT'))
111     elif str(anType) == 'Mjj':
112     mjj = True
113     systematics = eval(config.get('LimitGeneral','sys_Mjj'))
114     sys_cut_suffix=eval(config.get('LimitGeneral','sys_cut_suffix'))
115    
116    
117 nmohr 1.24
118 peller 1.23
119 peller 1.36 # -------------------- generate the Workspace with all Histograms: ----------------------------------------------------------------------
120 peller 1.11
121     WS = ROOT.RooWorkspace('%s'%Datacradbin,'%s'%Datacradbin) #Zee
122 peller 1.1 print 'WS initialized'
123 bortigno 1.7 disc= ROOT.RooRealVar(name,name,xMin,xMax)
124 peller 1.1 obs = ROOT.RooArgList(disc)
125     ROOT.gROOT.SetStyle("Plain")
126     datas = []
127     datatyps =[]
128     histos = []
129     typs = []
130 nmohr 1.37 hNames = []
131 nmohr 1.39 aNames = []
132 peller 1.1 statUps=[]
133     statDowns=[]
134 nmohr 1.15 if blind:
135 peller 1.29 printc('red','', 'I AM BLINDED!')
136 nmohr 1.15 counter=0
137 peller 1.29
138 peller 1.33 if weightF_sys:
139     #weightF_sys_function=config.get('Weights','weightF_sys')
140     weightF_sys_histos = []
141     weightF_sys_Ups = []
142     weightF_sys_Downs = []
143 nmohr 1.37 if addSample_sys:
144     sTyps = []
145     addSample_sys_histos = []
146     aSample_sys_Ups = []
147     aSample_sys_Downs = []
148 peller 1.33
149 peller 1.1 for job in info:
150 peller 1.28 if eval(job.active):
151 peller 1.29 if job.subsamples:
152     for subsample in range(0,len(job.subnames)):
153     if job.subnames[subsample] in BKGlist:
154 peller 1.45 print 'getting %s'%job.subnames[subsample]
155 nmohr 1.31 hTemp, typ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample)
156 peller 1.45 print hTemp.Integral()
157 peller 1.29 histos.append(hTemp)
158     typs.append(Group[job.subnames[subsample]])
159 nmohr 1.37 hNames.append(job.subnames[subsample])
160 peller 1.33 if weightF_sys:
161     hTempW, _ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample,'weightF_sys')
162     weightF_sys_histos.append(hTempW)
163 peller 1.29 if counter == 0:
164     hDummy = copy(hTemp)
165     else:
166     hDummy.Add(hTemp)
167     counter += 1
168    
169     elif job.subnames[subsample] == SIG:
170 peller 1.45 hNames.append(job.subnames[subsample])
171     print 'getting %s'%job.subnames[subsample]
172 nmohr 1.31 hTemp, typ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample)
173 peller 1.45 print hTemp.Integral()
174 peller 1.28 histos.append(hTemp)
175 peller 1.29 typs.append(Group[job.subnames[subsample]])
176 peller 1.33 if weightF_sys:
177     hTempW, _ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample,'weightF_sys')
178     weightF_sys_histos.append(hTempW)
179 nmohr 1.41 if addSample_sys and job.subnames[subsample] in addSample_sys.values():
180 nmohr 1.40 aNames.append(job.subnames[subsample])
181     hTempS, s_ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample)
182     addSample_sys_histos.append(hTempS)
183 peller 1.29
184     else:
185     if job.name in BKGlist:
186     #print job.getpath()
187 peller 1.45 print 'getting %s'%job.name
188 nmohr 1.31 hTemp, typ = getHistoFromTree(job,path,config,options,MC_rescale_factor)
189 peller 1.29 histos.append(hTemp)
190 peller 1.45 print hTemp.Integral()
191 peller 1.33 typs.append(Group[job.name])
192 peller 1.42 hNames.append(job.name)
193 peller 1.33 if weightF_sys:
194     hTempW, _ = getHistoFromTree(job,path,config,options,MC_rescale_factor,-1,'weightF_sys')
195     weightF_sys_histos.append(hTempW)
196    
197 peller 1.29 if counter == 0:
198     hDummy = copy(hTemp)
199     else:
200     hDummy.Add(hTemp)
201     counter += 1
202    
203     elif job.name == SIG:
204 peller 1.45 print 'getting %s'%job.name
205 nmohr 1.31 hTemp, typ = getHistoFromTree(job,path,config,options,MC_rescale_factor)
206 peller 1.28 histos.append(hTemp)
207 peller 1.45 print hTemp.Integral()
208 peller 1.33 typs.append(Group[job.name])
209 nmohr 1.37 hNames.append(job.name)
210 peller 1.33 if weightF_sys:
211     hTempW, _ = getHistoFromTree(job,path,config,options,MC_rescale_factor,-1,'weightF_sys')
212     weightF_sys_histos.append(hTempW)
213    
214 peller 1.29 elif job.name in data:
215     #print 'DATA'
216 peller 1.45 print 'getting %s'%job.name
217 nmohr 1.31 hTemp, typ = getHistoFromTree(job,path,config,options)
218 peller 1.29 datas.append(hTemp)
219 peller 1.45 print hTemp.Integral()
220 peller 1.29 datatyps.append(typ)
221 nmohr 1.40
222 peller 1.43 if addSample_sys and job.name in addSample_sys.values():
223 peller 1.42 aNames.append(job.name)
224     hTempS, s_ = getHistoFromTree(job,path,config,options,MC_rescale_factor)
225 nmohr 1.40 addSample_sys_histos.append(hTempS)
226 peller 1.1
227     MC_integral=0
228     MC_entries=0
229     for histo in histos:
230     MC_integral+=histo.Integral()
231 peller 1.45 print 'histo integral %s'%histo.Integral()
232 peller 1.36 printc('green','', 'MC integral = %s'%MC_integral)
233    
234 nmohr 1.37 def getAlternativeShapes(histos,altHistos,hNames,aNames,addSample_sys):
235     theHistosUp = copy(histos)
236     theHistosDown = copy(histos)
237     for name in addSample_sys.keys():
238     print name
239     hVar = altHistos[aNames.index(addSample_sys[name])].Clone()
240     hNom = histos[hNames.index(name)].Clone()
241     hAlt = hNom.Clone()
242     hNom.Add(hVar,-1.)
243     hAlt.Add(hNom)
244     for bin in range(0,nBins):
245     if hAlt.GetBinContent(bin) < 0.: hAlt.SetBinContent(bin,0.)
246     theHistosUp[hNames.index(name)] = hVar.Clone()
247     theHistosDown[hNames.index(name)] = hAlt.Clone()
248     return theHistosUp, theHistosDown
249    
250    
251    
252    
253 peller 1.13 #order and add together
254 peller 1.33 typs2=copy(typs)
255 nmohr 1.37 typs3=copy(typs)
256     typs4=copy(typs)
257 nmohr 1.40 if addSample_sys:
258     aSampleUp, aSampleDown = getAlternativeShapes(histos,addSample_sys_histos,hNames,aNames,addSample_sys)
259 peller 1.13 histos, typs = orderandadd(histos,typs,setup)
260 nmohr 1.37 if weightF_sys:
261     weightF_sys_histos,_=orderandadd(weightF_sys_histos,typs2,setup)
262     if addSample_sys:
263     aSampleUp,aNames=orderandadd(aSampleUp,typs3,setup)
264     aSampleDown,aNames=orderandadd(aSampleDown,typs4,setup)
265 peller 1.1
266     for i in range(0,len(histos)):
267 peller 1.29 newname=Dict[typs[i]]
268     histos[i].SetName(newname)
269 peller 1.5 #histos[i].SetDirectory(outfile)
270     outfile.cd()
271 peller 1.1 histos[i].Write()
272 nmohr 1.22 errorsum=0
273     total=0
274     for j in range(histos[i].GetNbinsX()+1):
275     errorsum=errorsum+(histos[i].GetBinError(j))**2
276     errorsum=sqrt(errorsum)
277     total=histos[i].Integral()
278 peller 1.35
279 peller 1.36 if binstat: #treating statistics in single bins
280 peller 1.35 for bin in range(0,nBins):
281     statUps.append(histos[i].Clone())
282     statDowns.append(histos[i].Clone())
283     statUps[i*nBins+bin].SetName('%sCMS_vhbb_stats_%s_%s_%sUp'%(newname,newname,bin,options[10]))
284     statDowns[i*nBins+bin].SetName('%sCMS_vhbb_stats_%s_%s_%sDown'%(newname,newname,bin,options[10]))
285     #shift up and down with statistical error
286     if rescaleSqrtN:
287     statUps[i*nBins+bin].SetBinContent(bin,statUps[i*nBins+bin].GetBinContent(bin)+statUps[i*nBins+bin].GetBinError(bin)/total*errorsum)
288     statDowns[i*nBins+bin].SetBinContent(bin,statDowns[i*nBins+bin].GetBinContent(bin)-statDowns[i*nBins+bin].GetBinError(bin)/total*errorsum)
289     else:
290     statUps[i*nBins+bin].SetBinContent(bin,statUps[i*nBins+bin].GetBinContent(bin)+statUps[i*nBins+bin].GetBinError(bin))
291     statDowns[i*nBins+bin].SetBinContent(bin,statDowns[i*nBins+bin].GetBinContent(bin)-statDowns[i*nBins+bin].GetBinError(bin))
292     statUps[i*nBins+bin].Write()
293     statDowns[i*nBins+bin].Write()
294     histPdf = ROOT.RooDataHist(newname,newname,obs,histos[i])
295     #UP stats of MCs
296     RooStatsUp = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%s_%sUp'%(newname,newname,bin,options[10]),'%sCMS_vhbb_stats_%s_%s_%sUp'%(newname,newname,bin,options[10]),obs, statUps[i*nBins+bin])
297     #DOWN stats of MCs
298     RooStatsDown = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%s_%sDown'%(newname,newname,bin,options[10]),'%sCMS_vhbb_stats_%s_%s_%sDown'%(newname,newname,bin,options[10]),obs, statDowns[i*nBins+bin])
299     getattr(WS,'import')(histPdf)
300     getattr(WS,'import')(RooStatsUp)
301     getattr(WS,'import')(RooStatsDown)
302 peller 1.13
303 peller 1.35 else:
304     statUps.append(histos[i].Clone())
305     statDowns.append(histos[i].Clone())
306     statUps[i].SetName('%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]))
307     statDowns[i].SetName('%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]))
308     #shift up and down with statistical error
309     for j in range(histos[i].GetNbinsX()+1):
310     if rescaleSqrtN:
311     statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j)/total*errorsum)
312     statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j)/total*errorsum)
313     else:
314     statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
315     statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
316 peller 1.47 #if statDowns[i].GetBinError(j)<0.: statDowns[i].SetBinContent(j,0.)
317 peller 1.35 statUps[i].Write()
318     statDowns[i].Write()
319     histPdf = ROOT.RooDataHist(newname,newname,obs,histos[i])
320     #UP stats of MCs
321     RooStatsUp = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]),'%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]),obs, statUps[i])
322     #DOWN stats of MCs
323     RooStatsDown = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]),'%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]),obs, statDowns[i])
324     getattr(WS,'import')(histPdf)
325     getattr(WS,'import')(RooStatsUp)
326     getattr(WS,'import')(RooStatsDown)
327    
328     #And now WeightF sys
329 peller 1.33 if weightF_sys:
330 peller 1.35 weightF_sys_Downs.append(weightF_sys_histos[i].Clone())
331     weightF_sys_Ups.append(weightF_sys_histos[i].Clone())
332     weightF_sys_Downs[i].SetName('%sCMS_vhbb_weightF_%sDown'%(newname,options[10]))
333     weightF_sys_Ups[i].SetName('%sCMS_vhbb_weightF_%sUp'%(newname,options[10]))
334     for j in range(histos[i].GetNbinsX()+1):
335     weightF_sys_Ups[i].SetBinContent(j,2*histos[i].GetBinContent(j)-weightF_sys_Downs[i].GetBinContent(j))
336 peller 1.33 weightF_sys_Ups[i].Write()
337 peller 1.35 weightF_sys_Downs[i].Write()
338 peller 1.33 RooWeightFUp = ROOT.RooDataHist('%sCMS_vhbb_weightF_%sUp'%(newname,options[10]),'%sCMS_vhbb_weightF_%s_%sUp'%(newname,newname,options[10]),obs, weightF_sys_Ups[i])
339     RooWeightFDown = ROOT.RooDataHist('%sCMS_vhbb_weightF_%sDown'%(newname,options[10]),'%sCMS_vhbb_weightF_%s_%sDown'%(newname,newname,options[10]),obs, weightF_sys_Downs[i])
340     getattr(WS,'import')(RooWeightFUp)
341     getattr(WS,'import')(RooWeightFDown)
342 nmohr 1.37 #And now Additional sample sys
343     if addSample_sys:
344     aSample_sys_Downs.append(aSampleUp[i].Clone())
345     aSample_sys_Ups.append(aSampleDown[i].Clone())
346     aSample_sys_Downs[i].SetName('%sCMS_vhbb_model_%sDown'%(newname,newname))
347     aSample_sys_Ups[i].SetName('%sCMS_vhbb_model_%sUp'%(newname,newname))
348     aSample_sys_Ups[i].Write()
349     aSample_sys_Downs[i].Write()
350     RooSampleUp = ROOT.RooDataHist('%sCMS_vhbb_model_%sUp'%(newname,newname),'%sCMS_vhbb_model_%sUp'%(newname,newname),obs, aSample_sys_Ups[i])
351     RooSampleDown = ROOT.RooDataHist('%sCMS_vhbb_model_%sDown'%(newname,newname),'%sCMS_vhbb_model_%sDown'%(newname,newname),obs, aSample_sys_Downs[i])
352     getattr(WS,'import')(RooSampleUp)
353     getattr(WS,'import')(RooSampleDown)
354 peller 1.1
355    
356     #HISTOGRAMM of DATA
357     d1 = ROOT.TH1F('d1','d1',nBins,xMin,xMax)
358     for i in range(0,len(datas)):
359     d1.Add(datas[i],1)
360 peller 1.13 printc('green','','\nDATA integral = %s\n'%d1.Integral())
361 peller 1.1 flow = d1.GetEntries()-d1.Integral()
362     if flow > 0:
363 peller 1.13 printc('red','','U/O flow: %s'%flow)
364 peller 1.29 d1.SetName(Dict['Data'])
365 peller 1.5 outfile.cd()
366 peller 1.47 #d1.Write()
367    
368 peller 1.36
369 nmohr 1.15 if blind:
370 peller 1.47 print 'toy data integral: %s'%hDummy.Integral()
371 peller 1.29 hDummy.SetName(Dict['Data'])
372     histPdf = ROOT.RooDataHist(Dict['Data'],Dict['Data'],obs,hDummy)
373 nmohr 1.25 #rooDummy = ROOT.RooDataHist('data_obs','data_obs',obs,hDummy)
374     #toy = ROOT.RooHistPdf('data_obs','data_obs',ROOT.RooArgSet(obs),rooDummy)
375     #rooDataSet = toy.generate(ROOT.RooArgSet(obs),int(d1.Integral()))
376     #histPdf = ROOT.RooDataHist('data_obs','data_obs',ROOT.RooArgSet(obs),rooDataSet.reduce(ROOT.RooArgSet(obs)))
377 peller 1.47 hDummy.Write()
378 nmohr 1.15 else:
379 peller 1.29 histPdf = ROOT.RooDataHist(Dict['Data'],Dict['Data'],obs,d1)
380 peller 1.47 d1.Write()
381 peller 1.1 #ROOT.RooAbsData.plotOn(histPdf,frame)
382     getattr(WS,'import')(histPdf)
383    
384     #SYSTEMATICS:
385     UD = ['Up','Down']
386     systhistosarray=[]
387 peller 1.13 Coco=0 #iterates over (all systematics) * (up,down)
388 nmohr 1.22 nominalShape = options[0]
389 bortigno 1.7
390     for sys in systematics:
391 peller 1.36 for Q in UD: # Q = 'Up' and 'Down'
392 peller 1.32 #options[7] ist der CutString name
393     new_cut=sys_cut_suffix[sys]
394     new_options = copy(options)
395     if not new_cut == 'nominal':
396 peller 1.34 old_str,new_str=new_cut.split('>')
397     new_options[7]=[options[7],old_str,new_str.replace('?',Q)]
398 peller 1.1 ff=options[0].split('.')
399 bortigno 1.7 if bdt == True:
400 peller 1.13 ff[1]='%s_%s'%(sys,Q.lower())
401 peller 1.32 new_options[0]=nominalShape.replace('.nominal','.%s_%s'%(sys,Q.lower()))
402 bortigno 1.7 elif mjj == True:
403     ff[0]='H_%s'%(sys)
404 peller 1.13 ff[1]='mass_%s'%(Q.lower())
405 peller 1.32 new_options[0]='.'.join(ff)
406 peller 1.1
407 peller 1.8 print '\n'
408 peller 1.13 printc('blue','','\t--> doing systematic %s %s'%(sys,Q.lower()))
409 peller 1.1
410     systhistosarray.append([])
411     typsX = []
412    
413     for job in info:
414 peller 1.29 if eval(job.active):
415     if job.subsamples:
416     for subsample in range(0,len(job.subnames)):
417     if job.subnames[subsample] in BKGlist:
418 peller 1.32 hTemp, typ = getHistoFromTree(job,path,config,new_options,MC_rescale_factor,subsample)
419 peller 1.29 systhistosarray[Coco].append(hTemp)
420     typsX.append(Group[job.subnames[subsample]])
421     elif job.subnames[subsample] == SIG:
422 peller 1.32 hTemp, typ = getHistoFromTree(job,path,config,new_options,MC_rescale_factor,subsample)
423 peller 1.29 systhistosarray[Coco].append(hTemp)
424     typsX.append(Group[job.subnames[subsample]])
425    
426     else:
427     if job.name in BKGlist:
428 peller 1.32 hTemp, typ = getHistoFromTree(job,path,config,new_options,MC_rescale_factor)
429 peller 1.29 systhistosarray[Coco].append(hTemp)
430     typsX.append(Group[job.name])
431     elif job.name == SIG:
432 peller 1.32 hTemp, typ = getHistoFromTree(job,path,config,new_options,MC_rescale_factor)
433 peller 1.29 systhistosarray[Coco].append(hTemp)
434     typsX.append(Group[job.name])
435    
436 peller 1.1 MC_integral=0
437     for histoX in systhistosarray[Coco]:
438     MC_integral+=histoX.Integral()
439 peller 1.13 printc('green','', 'MC integral = %s'%MC_integral)
440     systhistosarray[Coco], typsX = orderandadd(systhistosarray[Coco],typsX,setup)
441 peller 1.29
442 peller 1.36 if scaling: #rescaling after the sys has been propagated through the BDT with a scaling
443 peller 1.13 for i in range(0,len(systhistosarray[Coco])):
444     for bin in range(0,histos[i].GetSize()):
445     A=systhistosarray[Coco][i].GetBinContent(bin)
446     B=histos[i].GetBinContent(bin)
447     systhistosarray[Coco][i].SetBinContent(bin,B+((A-B)/4.))
448 nmohr 1.24 # finaly lpop over histos
449 peller 1.1 for i in range(0,len(systhistosarray[Coco])):
450 peller 1.29 systhistosarray[Coco][i].SetName('%s%s%s'%(Dict[typs[i]],systematicsnaming[sys],Q))
451 peller 1.5 outfile.cd()
452 peller 1.13 systhistosarray[Coco][i].Write()
453 peller 1.29 histPdf = ROOT.RooDataHist('%s%s%s'%(Dict[typs[i]],systematicsnaming[sys],Q),'%s%s%s'%(Dict[typs[i]],systematicsnaming[sys],Q),obs,systhistosarray[Coco][i])
454 peller 1.1 getattr(WS,'import')(histPdf)
455     Coco+=1
456     WS.writeToFile(outpath+'vhbb_WS_'+ROOToutname+'.root')
457    
458 peller 1.36
459    
460     # -------------------- write DATAcard: ----------------------------------------------------------------------
461 peller 1.29 columns=len(setup)
462    
463 nmohr 1.18 if '8TeV' in options[10]:
464 peller 1.42 pier = open(vhbbpath+'/python/pier8TeV.txt','r')
465 nmohr 1.18 else:
466 peller 1.42 pier = open(vhbbpath+'/python/pier.txt','r')
467 peller 1.1 scalefactors=pier.readlines()
468     pier.close()
469     f = open(outpath+'vhbb_DC_'+ROOToutname+'.txt','w')
470     f.write('imax\t1\tnumber of channels\n')
471 peller 1.29 f.write('jmax\t%s\tnumber of backgrounds (\'*\' = automatic)\n'%(columns-1))
472 peller 1.1 f.write('kmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
473 peller 1.11 if bdt==True:
474     f.write('shapes * * vhbb_WS_%s.root $CHANNEL:$PROCESS $CHANNEL:$PROCESS$SYSTEMATIC\n\n'%ROOToutname)
475     else:
476     f.write('shapes * * vhbb_TH_%s.root $PROCESS $PROCESS$SYSTEMATIC\n\n'%ROOToutname)
477     f.write('bin\t%s\n\n'%Datacradbin)
478 nmohr 1.25 if blind:
479     f.write('observation\t%s\n\n'%(hDummy.Integral()))
480     else:
481     f.write('observation\t%s\n\n'%(int(d1.Integral())))
482 peller 1.14
483 peller 1.29 f.write('bin')
484     for c in range(0,columns): f.write('\t%s'%Datacradbin)
485     f.write('\n')
486    
487     f.write('process')
488     for c in setup: f.write('\t%s'%Dict[c])
489     f.write('\n')
490    
491     f.write('process')
492     for c in range(0,columns): f.write('\t%s'%c)
493     f.write('\n')
494    
495     f.write('rate')
496     for c in range(0,columns): f.write('\t%s'%histos[c].Integral())
497     f.write('\n')
498    
499     InUse=eval(config.get('Datacard','InUse'))
500     #Parse from config
501     for item in InUse:
502     f.write(item)
503     what=eval(config.get('Datacard',item))
504     f.write('\t%s'%what['type'])
505     for c in setup:
506     if c in what:
507     if item == 'CMS_eff_e' and 'Zmm' in options[10]: f.write('\t-')
508     elif item == 'CMS_eff_m' and 'Zee' in options[10]: f.write('\t-')
509     elif item == 'CMS_trigger_e' and 'Zmm' in options[10]: f.write('\t-')
510     elif item == 'CMS_trigger_m' and 'Zee' in options[10]: f.write('\t-')
511     else:
512     f.write('\t%s'%what[c])
513     else:
514     f.write('\t-')
515     f.write('\n')
516    
517     #Write shape stats and sys
518 peller 1.35 if binstat:
519     for c in setup:
520     for bin in range(0,nBins):
521     f.write('CMS_vhbb_stats_%s_%s_%s\tshape'%(Dict[c], bin, options[10]))
522     for it in range(0,columns):
523     if it == setup.index(c):
524     f.write('\t1.0')
525     else:
526     f.write('\t-')
527     f.write('\n')
528    
529     else:
530     for c in setup:
531     f.write('CMS_vhbb_stats_%s_%s\tshape'%(Dict[c], options[10]))
532     for it in range(0,columns):
533     if it == setup.index(c):
534     f.write('\t1.0')
535     else:
536     f.write('\t-')
537     f.write('\n')
538 peller 1.29
539 peller 1.33 if weightF_sys:
540     f.write('CMS_vhbb_weightF_%s\tshape'%(options[10]))
541     for it in range(0,columns): f.write('\t1.0')
542     f.write('\n')
543    
544 nmohr 1.37 if addSample_sys:
545     for newSample in addSample_sys.iterkeys():
546     for c in setup:
547 nmohr 1.39 if not c == Group[newSample]: continue
548 nmohr 1.37 f.write('CMS_vhbb_model_%s\tshape'%(Dict[c]))
549     for it in range(0,columns):
550     if it == setup.index(c):
551     f.write('\t1.0')
552     else:
553     f.write('\t-')
554     f.write('\n')
555 peller 1.33
556 peller 1.45 if scaling: sys_factor=1.0
557 peller 1.29 else: sys_factor=1.0
558     for sys in systematics:
559     f.write('%s\tshape'%systematicsnaming[sys])
560     for c in range(0,columns): f.write('\t%s'%sys_factor)
561     f.write('\n')
562 peller 1.1 f.close()
563 nmohr 1.31 outfile.Close()