ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/workspace_datacard.py
Revision: 1.93
Committed: Tue Apr 23 13:00:57 2013 UTC (12 years ago) by nmohr
Content type: text/x-python
Branch: MAIN
CVS Tags: lhcp_UnblindFix, hcp_Unblind, HEAD
Changes since 1.92: +19 -2 lines
Log Message:
Revive signal inject

File Contents

# User Rev Content
1 peller 1.8 #!/usr/bin/env python
2 peller 1.53 import os, sys, ROOT, warnings, pickle
3 nmohr 1.68 ROOT.gROOT.SetBatch(True)
4 peller 1.1 from array import array
5     from math import sqrt
6 nmohr 1.68 from copy import copy, deepcopy
7 peller 1.1 #suppres the EvalInstace conversion warning bug
8     warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
9 nmohr 1.31 from optparse import OptionParser
10 nmohr 1.80 from myutils import BetterConfigParser, Sample, progbar, printc, ParseInfo, Rebinner, HistoMaker
11 peller 1.1
12 nmohr 1.68 #--CONFIGURE---------------------------------------------------------------------
13 nmohr 1.31 argv = sys.argv
14     parser = OptionParser()
15 bortigno 1.79 parser.add_option("-V", "--variable", dest="variable", default="",
16 nmohr 1.31 help="variable for shape analysis")
17     parser.add_option("-C", "--config", dest="config", default=[], action="append",
18     help="configuration file")
19     (opts, args) = parser.parse_args(argv)
20 nmohr 1.20 config = BetterConfigParser()
21 nmohr 1.31 config.read(opts.config)
22 nmohr 1.68 var=opts.variable
23     #-------------------------------------------------------------------------------
24    
25     #--read variables from config---------------------------------------------------
26     # 7 or 8TeV Analysis
27 nmohr 1.31 anaTag = config.get("Analysis","tag")
28 nmohr 1.68 if not any([anaTag == '7TeV',anaTag == '8TeV']):
29     raise Exception("anaTag %s unknown. Specify 8TeV or 7TeV in the general config"%anaTag)
30     # Directories:
31 peller 1.1 Wdir=config.get('Directories','Wdir')
32 peller 1.42 vhbbpath=config.get('Directories','vhbbpath')
33 nmohr 1.46 samplesinfo=config.get('Directories','samplesinfo')
34 peller 1.67 path = config.get('Directories','dcSamples')
35 nmohr 1.68 outpath=config.get('Directories','limits')
36     try:
37     os.stat(outpath)
38     except:
39     os.mkdir(outpath)
40     # parse histogram config:
41 nmohr 1.80 treevar = config.get('dc:%s'%var,'var')
42     name = config.get('dc:%s'%var,'wsVarName')
43     title = name
44     nBins = int(config.get('dc:%s'%var,'range').split(',')[0])
45     xMin = float(config.get('dc:%s'%var,'range').split(',')[1])
46     xMax = float(config.get('dc:%s'%var,'range').split(',')[2])
47     ROOToutname = config.get('dc:%s'%var,'dcName')
48     RCut = config.get('dc:%s'%var,'cut')
49     signals = config.get('dc:%s'%var,'signal').split(' ')
50     datas = config.get('dc:%s'%var,'dcBin')
51     Datacardbin=config.get('dc:%s'%var,'dcBin')
52     anType = config.get('dc:%s'%var,'type')
53 peller 1.29 setup=eval(config.get('LimitGeneral','setup'))
54 nmohr 1.68 #Systematics:
55     if config.has_option('LimitGeneral','addSample_sys'):
56     addSample_sys = eval(config.get('LimitGeneral','addSample_sys'))
57     additionals = [addSample_sys[key] for key in addSample_sys]
58     else:
59     addSample_sys = None
60     additionals = []
61     #find out if BDT or MJJ:
62     bdt = False
63     mjj = False
64 bortigno 1.78 cr = False
65 nmohr 1.68 if str(anType) == 'BDT':
66     bdt = True
67     systematics = eval(config.get('LimitGeneral','sys_BDT'))
68     elif str(anType) == 'Mjj':
69     mjj = True
70     systematics = eval(config.get('LimitGeneral','sys_Mjj'))
71 bortigno 1.78 elif str(anType) == 'cr':
72     cr = True
73     systematics = eval(config.get('LimitGeneral','sys_cr'))
74    
75 nmohr 1.68 sys_cut_suffix=eval(config.get('LimitGeneral','sys_cut_suffix'))
76 peller 1.86 sys_cut_include=[]
77     if config.has_option('LimitGeneral','sys_cut_include'):
78     sys_cut_include=eval(config.get('LimitGeneral','sys_cut_include'))
79 nmohr 1.68 systematicsnaming = eval(config.get('LimitGeneral','systematicsnaming'))
80     sys_factor_dict = eval(config.get('LimitGeneral','sys_factor'))
81     sys_affecting = eval(config.get('LimitGeneral','sys_affecting'))
82     # weightF:
83     weightF = config.get('Weights','weightF')
84 nmohr 1.84 weightF_systematics = eval(config.get('LimitGeneral','weightF_sys'))
85 peller 1.77 # rescale stat shapes by sqrtN
86     rescaleSqrtN=eval(config.get('LimitGeneral','rescaleSqrtN'))
87 nmohr 1.68 # get nominal cutstring:
88     treecut = config.get('Cuts',RCut)
89     # Train flag: splitting of samples
90     TrainFlag = eval(config.get('Analysis','TrainFlag'))
91 peller 1.75 # toy data option:
92     toy=eval(config.get('LimitGeneral','toy'))
93 nmohr 1.68 # blind data option:
94     blind=eval(config.get('LimitGeneral','blind'))
95 nmohr 1.81 # additional blinding cut:
96     addBlindingCut = None
97     if config.has_option('LimitGeneral','addBlindingCut'):
98     addBlindingCut = config.get('LimitGeneral','addBlindingCut')
99 peller 1.89 print 'adding add. blinding cut'
100 peller 1.83 #change nominal shapes by syst
101     change_shapes = None
102     if config.has_option('LimitGeneral','change_shapes'):
103     change_shapes = eval(config.get('LimitGeneral','change_shapes'))
104     print 'changing the shapes'
105 bortigno 1.78 #on control region cr never blind. Overwrite whatever is in the config
106     if str(anType) == 'cr':
107     if blind:
108     print '@WARNING: Changing blind to false since you are running for control region.'
109     blind = False
110 nmohr 1.68 if blind:
111 bortigno 1.78 printc('red','', 'I AM BLINDED!')
112 nmohr 1.68 #get List of backgrounds in use:
113     backgrounds = eval(config.get('LimitGeneral','BKG'))
114     #Groups for adding samples together
115 peller 1.75 GroupDict = eval(config.get('LimitGeneral','Group'))
116 nmohr 1.68 #naming for DC
117     Dict= eval(config.get('LimitGeneral','Dict'))
118     #treating statistics bin-by-bin:
119     binstat = eval(config.get('LimitGeneral','binstat'))
120     # Use the rebinning:
121     rebin_active=eval(config.get('LimitGeneral','rebin_active'))
122 bortigno 1.78 if str(anType) == 'cr':
123     if rebin_active:
124     print '@WARNING: Changing rebin_active to false since you are running for control region.'
125     rebin_active = False
126 peller 1.77 # ignore stat shapes
127     ignore_stats = eval(config.get('LimitGeneral','ignore_stats'))
128 peller 1.72 #max_rel = float(config.get('LimitGeneral','rebin_max_rel'))
129 nmohr 1.68 signal_inject=config.get('LimitGeneral','signal_inject')
130     # add signal as background
131     add_signal_as_bkg=config.get('LimitGeneral','add_signal_as_bkg')
132     if not add_signal_as_bkg == 'None':
133     setup.append(add_signal_as_bkg)
134     #----------------------------------------------------------------------------
135 peller 1.65
136 nmohr 1.68 #--Setup--------------------------------------------------------------------
137     #Assign Pt region for sys factors
138 peller 1.58 if 'HighPtLooseBTag' in ROOToutname:
139     pt_region = 'HighPtLooseBTag'
140 nmohr 1.90 elif 'HighPt' in ROOToutname or 'highPt' in ROOToutname:
141 peller 1.58 pt_region = 'HighPt'
142 nmohr 1.90 elif 'MedPt' in ROOToutname:
143     pt_region = 'MedPt'
144 peller 1.60 elif 'LowPt' in ROOToutname or 'lowPt' in ROOToutname:
145 peller 1.58 pt_region = 'LowPt'
146 nmohr 1.68 elif 'ATLAS' in ROOToutname:
147     pt_region = 'HighPt'
148 nmohr 1.80 elif 'MJJ' in ROOToutname:
149 nmohr 1.68 pt_region = 'HighPt'
150     else:
151 peller 1.58 print "Unknown Pt region"
152     sys.exit("Unknown Pt region")
153 nmohr 1.68 # Set rescale factor of 2 in case of TrainFalg
154 peller 1.44 if TrainFlag:
155     MC_rescale_factor=2.
156     print 'I RESCALE BY 2.0'
157     else: MC_rescale_factor = 1.
158 nmohr 1.68 #systematics up/down
159     UD = ['Up','Down']
160     #Parse samples configuration
161     info = ParseInfo(samplesinfo,path)
162     # get all the treeCut sets
163     # create different sample Lists
164     all_samples = info.get_samples(signals+backgrounds+additionals)
165 nmohr 1.93
166 nmohr 1.68 signal_samples = info.get_samples(signals)
167     background_samples = info.get_samples(backgrounds)
168 nmohr 1.80 data_sample_names = config.get('dc:%s'%var,'data').split(' ')
169 peller 1.77 data_samples = info.get_samples(data_sample_names)
170 nmohr 1.68 #-------------------------------------------------------------------------------------------------
171    
172     optionsList=[]
173    
174 nmohr 1.82 def appendList(): optionsList.append({'cut':copy(_cut),'var':copy(_treevar),'name':copy(_name),'nBins':nBins,'xMin':xMin,'xMax':xMax,'weight':copy(_weight),'blind':blind})
175 nmohr 1.68
176     #nominal
177     _cut = treecut
178     _treevar = treevar
179     _name = title
180     _weight = weightF
181     appendList()
182    
183     #the 4 sys
184     for syst in systematics:
185     for Q in UD:
186     #default:
187     _cut = treecut
188     _name = title
189     _weight = weightF
190     #replace cut string
191     new_cut=sys_cut_suffix[syst]
192     if not new_cut == 'nominal':
193     old_str,new_str=new_cut.split('>')
194     _cut = treecut.replace(old_str,new_str.replace('?',Q))
195     _name = title
196     _weight = weightF
197     #replace tree variable
198     if bdt == True:
199 peller 1.77 #ff[1]='%s_%s'%(sys,Q.lower())
200     _treevar = treevar.replace('.nominal','.%s_%s'%(syst,Q.lower()))
201     print _treevar
202 nmohr 1.68 elif mjj == True:
203 nmohr 1.80 if syst == 'JER' or syst == 'JES':
204     _treevar = 'H_%s.mass_%s'%(syst,Q.lower())
205 nmohr 1.68 else:
206     _treevar = treevar
207 bortigno 1.78 elif cr == True:
208 nmohr 1.80 if syst == 'beff' or syst == 'bmis' or syst == 'beff1':
209     _treevar = treevar.replace(old_str,new_str.replace('?',Q))
210 bortigno 1.78 else:
211     _treevar = treevar
212 nmohr 1.68 #append
213     appendList()
214    
215     #UEPS
216 nmohr 1.84 for weightF_sys in weightF_systematics:
217     for _weight in [config.get('Weights','%s_UP' %(weightF_sys)),config.get('Weights','%s_DOWN' %(weightF_sys))]:
218 nmohr 1.68 _cut = treecut
219     _treevar = treevar
220     _name = title
221     appendList()
222    
223     #for option in optionsList:
224     # print option['cut']
225    
226 peller 1.75 mc_hMaker = HistoMaker(all_samples,path,config,optionsList,GroupDict)
227 nmohr 1.68 data_hMaker = HistoMaker(data_samples,path,config,[optionsList[0]])
228     #Calculate lumi
229     lumi = 0.
230     nData = 0
231     for job in data_samples:
232     nData += 1
233     lumi += float(job.lumi)
234 peller 1.57
235 nmohr 1.68 if nData > 1:
236     lumi = lumi/float(nData)
237 nmohr 1.24
238 nmohr 1.68 mc_hMaker.lumi = lumi
239     data_hMaker.lumi = lumi
240 peller 1.62
241 nmohr 1.81 if addBlindingCut:
242     for i in range(len(mc_hMaker.optionsList)):
243     mc_hMaker.optionsList[i]['cut'] += ' & %s' %addBlindingCut
244     for i in range(len(data_hMaker.optionsList)):
245     data_hMaker.optionsList[i]['cut'] += ' & %s' %addBlindingCut
246 peller 1.69
247 nmohr 1.82
248 peller 1.72 if rebin_active:
249     mc_hMaker.calc_rebin(background_samples)
250     #transfer rebinning info to data maker
251 peller 1.74 data_hMaker.norebin_nBins = copy(mc_hMaker.norebin_nBins)
252     data_hMaker.rebin_nBins = copy(mc_hMaker.rebin_nBins)
253 peller 1.77 data_hMaker.nBins = copy(mc_hMaker.nBins)
254     data_hMaker._rebin = copy(mc_hMaker._rebin)
255 peller 1.74 data_hMaker.mybinning = deepcopy(mc_hMaker.mybinning)
256 peller 1.69
257 nmohr 1.68 all_histos = {}
258     data_histos = {}
259 peller 1.53
260 peller 1.75 print '\n\t...fetching histos...'
261    
262 nmohr 1.68 for job in all_samples:
263 peller 1.75 print '\t- %s'%job
264 nmohr 1.88 if not GroupDict[job.name] in sys_cut_include:
265 peller 1.86 # manual overwrite
266 peller 1.89 if addBlindingCut:
267     all_histos[job.name] = mc_hMaker.get_histos_from_tree(job,treecut+'& %s'%addBlindingCut)
268     else:
269     all_histos[job.name] = mc_hMaker.get_histos_from_tree(job,treecut)
270 peller 1.86 else:
271     all_histos[job.name] = mc_hMaker.get_histos_from_tree(job)
272 peller 1.53
273 nmohr 1.68 for job in data_samples:
274 peller 1.75 print '\t- %s'%job
275 nmohr 1.68 data_histos[job.name] = data_hMaker.get_histos_from_tree(job)[0]['DATA']
276 peller 1.23
277 peller 1.75 print '\t> done <\n'
278    
279     i=0
280     for job in background_samples:
281     print job.name
282     htree = all_histos[job.name][0].values()[0]
283     if not i:
284     hDummy = copy(htree)
285     else:
286     hDummy.Add(htree,1)
287     del htree
288     i+=1
289    
290 nmohr 1.93 if signal_inject:
291     signal_inject = info.get_samples([signal_inject])
292     sig_hMaker = HistoMaker(signal_inject,path,config,optionsList,GroupDict)
293     sig_hMaker.lumi = lumi
294     if rebin_active:
295     sig_hMaker.norebin_nBins = copy(mc_hMaker.norebin_nBins)
296     sig_hMaker.rebin_nBins = copy(mc_hMaker.rebin_nBins)
297     sig_hMaker.nBins = copy(mc_hMaker.nBins)
298     sig_hMaker._rebin = copy(mc_hMaker._rebin)
299     sig_hMaker.mybinning = deepcopy(mc_hMaker.mybinning)
300    
301     for job in signal_inject:
302     htree = sig_hMaker.get_histos_from_tree(job)
303     hDummy.Add(htree[0].values()[0],1)
304     del htree
305    
306 nmohr 1.68 nData = 0
307     for job in data_histos:
308     if nData == 0:
309     theData = data_histos[job]
310     else:
311     theData.Add(data_histos[i])
312 peller 1.11
313 nmohr 1.68 #-- Write Files-----------------------------------------------------------------------------------
314     # generate the TH outfile:
315     outfile = ROOT.TFile(outpath+'vhbb_TH_'+ROOToutname+'.root', 'RECREATE')
316     outfile.mkdir(Datacardbin,Datacardbin)
317     outfile.cd(Datacardbin)
318     # generate the Workspace:
319 peller 1.62 WS = ROOT.RooWorkspace('%s'%Datacardbin,'%s'%Datacardbin) #Zee
320 peller 1.1 print 'WS initialized'
321 bortigno 1.7 disc= ROOT.RooRealVar(name,name,xMin,xMax)
322 peller 1.1 obs = ROOT.RooArgList(disc)
323 nmohr 1.68 #
324 peller 1.1 ROOT.gROOT.SetStyle("Plain")
325 peller 1.53
326 nmohr 1.68 #order and add all together
327     final_histos = {}
328 peller 1.57
329 nmohr 1.68 print '\n\t--> Ordering and Adding Histos\n'
330 peller 1.53
331 nmohr 1.68 #NOMINAL:
332     final_histos['nominal'] = HistoMaker.orderandadd([all_histos['%s'%job][0] for job in all_samples],setup)
333 peller 1.53
334 nmohr 1.68 #SYSTEMATICS:
335     ind = 1
336     for syst in systematics:
337     for Q in UD:
338 nmohr 1.71 final_histos['%s_%s'%(systematicsnaming[syst],Q)] = HistoMaker.orderandadd([all_histos[job.name][ind] for job in all_samples],setup)
339 nmohr 1.68 ind+=1
340 nmohr 1.84 for weightF_sys in weightF_systematics:
341 nmohr 1.68 for Q in UD:
342 nmohr 1.84 final_histos['%s_%s'%(systematicsnaming[weightF_sys],Q)]= HistoMaker.orderandadd([all_histos[job.name][ind] for job in all_samples],setup)
343 nmohr 1.68 ind+=1
344    
345 peller 1.83 if change_shapes:
346     for key in change_shapes:
347     syst,val=change_shapes[key].split('*')
348 peller 1.85 final_histos[syst][key].Scale(float(val))
349     print 'scaled %s times %s val'%(syst,val)
350 peller 1.83
351    
352 nmohr 1.71 def get_alternate_shape(hNominal,hAlternate):
353     hVar = hAlternate.Clone()
354     hNom = hNominal.Clone()
355     hAlt = hNom.Clone()
356     hNom.Add(hVar,-1.)
357     hAlt.Add(hNom)
358     for bin in range(0,hNominal.GetNbinsX()+1):
359     if hAlt.GetBinContent(bin) < 0.: hAlt.SetBinContent(bin,0.)
360     return hVar,hAlt
361    
362     def get_alternate_shapes(all_histos,asample_dict,all_samples):
363     alternate_shapes_up = []
364     alternate_shapes_down = []
365     for job in all_samples:
366     nominal = all_histos[job.name][0]
367     if job.name in asample_dict:
368     alternate = copy(all_histos[asample_dict[job.name]][0])
369     hUp, hDown = get_alternate_shape(nominal[nominal.keys()[0]],alternate[alternate.keys()[0]])
370     alternate_shapes_up.append({nominal.keys()[0]:hUp})
371     alternate_shapes_down.append({nominal.keys()[0]:hDown})
372     else:
373 peller 1.74 newh=nominal[nominal.keys()[0]].Clone('%s_%s_Up'%(nominal[nominal.keys()[0]].GetName(),'model'))
374     alternate_shapes_up.append({nominal.keys()[0]:nominal[nominal.keys()[0]].Clone()})
375     alternate_shapes_down.append({nominal.keys()[0]:nominal[nominal.keys()[0]].Clone()})
376 nmohr 1.71 return alternate_shapes_up, alternate_shapes_down
377    
378     if addSample_sys:
379     aUp, aDown = get_alternate_shapes(all_histos,addSample_sys,all_samples)
380     final_histos['%s_Up'%(systematicsnaming['model'])]= HistoMaker.orderandadd(aUp,setup)
381 peller 1.74 del aUp
382 nmohr 1.71 final_histos['%s_Down'%(systematicsnaming['model'])]= HistoMaker.orderandadd(aDown,setup)
383    
384 peller 1.77
385     if not ignore_stats:
386     #make statistical shapes:
387 peller 1.87 if not binstat:
388 peller 1.77 for Q in UD:
389 peller 1.87 final_histos['%s_%s'%(systematicsnaming['stats'],Q)] = {}
390     for job,hist in final_histos['nominal'].items():
391     errorsum=0
392 peller 1.77 for j in range(hist.GetNbinsX()+1):
393 peller 1.87 errorsum=errorsum+(hist.GetBinError(j))**2
394     errorsum=sqrt(errorsum)
395     total=hist.Integral()
396     for Q in UD:
397     final_histos['%s_%s'%(systematicsnaming['stats'],Q)][job] = hist.Clone()
398     for j in range(hist.GetNbinsX()+1):
399     if Q == 'Up':
400     if rescaleSqrtN and total:
401     final_histos['%s_%s'%(systematicsnaming['stats'],Q)][job].SetBinContent(j,max(0,hist.GetBinContent(j)+hist.GetBinError(j)/total*errorsum))
402     else:
403     final_histos['%s_%s'%(systematicsnaming['stats'],Q)][job].SetBinContent(j,max(0,hist.GetBinContent(j)+hist.GetBinError(j)))
404     if Q == 'Down':
405     if rescaleSqrtN and total:
406     final_histos['%s_%s'%(systematicsnaming['stats'],Q)][job].SetBinContent(j,max(0,hist.GetBinContent(j)-hist.GetBinError(j)/total*errorsum))
407     else:
408     final_histos['%s_%s'%(systematicsnaming['stats'],Q)][job].SetBinContent(j,max(0,hist.GetBinContent(j)-hist.GetBinError(j)))
409     else:
410 nmohr 1.91 binsBelowThreshold = {}
411 peller 1.87 for bin in range(0,nBins):
412     for Q in UD:
413     final_histos['%s_bin%s_%s'%(systematicsnaming['stats'],bin,Q)] = {}
414     for job,hist in final_histos['nominal'].items():
415 nmohr 1.91 binsBelowThreshold[job] = []
416     if hist.GetBinContent(bin) > 0.:
417     if hist.GetBinError(bin)/sqrt(hist.GetBinContent(bin)) > 0.5 and hist.GetBinContent(bin) >= 1.:
418     binsBelowThreshold[job].append(bin)
419     elif hist.GetBinError(bin)/(hist.GetBinContent(bin)) > 0.5 and hist.GetBinContent(bin) < 1.:
420     binsBelowThreshold[job].append(bin)
421 peller 1.87 for Q in UD:
422     final_histos['%s_bin%s_%s'%(systematicsnaming['stats'],bin,Q)][job] = hist.Clone()
423     if Q == 'Up':
424     final_histos['%s_bin%s_%s'%(systematicsnaming['stats'],bin,Q)][job].SetBinContent(bin,max(0,hist.GetBinContent(bin)+hist.GetBinError(bin)))
425     if Q == 'Down':
426     final_histos['%s_bin%s_%s'%(systematicsnaming['stats'],bin,Q)][job].SetBinContent(bin,max(0,hist.GetBinContent(bin)-hist.GetBinError(bin)))
427 peller 1.83
428    
429 nmohr 1.68 #write shapes in WS:
430     for key in final_histos:
431     for job, hist in final_histos[key].items():
432     if 'nominal' == key:
433 nmohr 1.71 hist.SetName('%s'%(Dict[job]))
434 peller 1.70 hist.Write()
435 nmohr 1.71 rooDataHist = ROOT.RooDataHist('%s' %(Dict[job]),'%s'%(Dict[job]),obs, hist)
436 nmohr 1.68 getattr(WS,'import')(rooDataHist)
437     for Q in UD:
438     if Q in key:
439     theSyst = key.replace('_%s'%Q,'')
440     else:
441     continue
442 nmohr 1.71 if systematicsnaming['stats'] in key:
443     nameSyst = '%s_%s_%s' %(theSyst,Dict[job],Datacardbin)
444     elif systematicsnaming['model'] in key:
445     nameSyst = '%s_%s' %(theSyst,Dict[job])
446 nmohr 1.68 else:
447     nameSyst = theSyst
448 nmohr 1.71 hist.SetName('%s%s%s' %(Dict[job],nameSyst,Q))
449 peller 1.70 hist.Write()
450 nmohr 1.71 rooDataHist = ROOT.RooDataHist('%s%s%s' %(Dict[job],nameSyst,Q),'%s%s%s'%(Dict[job],nameSyst,Q),obs, hist)
451 nmohr 1.68 getattr(WS,'import')(rooDataHist)
452 peller 1.48
453 nmohr 1.93 if toy or signal_inject:
454 nmohr 1.82 hDummy.SetName('data_obs')
455     hDummy.Write()
456     rooDataHist = ROOT.RooDataHist('data_obs','data_obs',obs, hDummy)
457 peller 1.75 else:
458 nmohr 1.82 theData.SetName('data_obs')
459     theData.Write()
460     rooDataHist = ROOT.RooDataHist('data_obs','data_obs',obs, theData)
461 peller 1.75
462 nmohr 1.68 getattr(WS,'import')(rooDataHist)
463 peller 1.13
464 nmohr 1.68 WS.writeToFile(outpath+'vhbb_WS_'+ROOToutname+'.root')
465 peller 1.47
466 nmohr 1.68 # now we have a Dict final_histos with sets of all grouped MCs for all systematics:
467     # nominal, ($SYS_Up/Down)*4, weightF_sys_Up/Down, stats_Up/Down
468 peller 1.36
469 peller 1.69 print '\n\t >>> PRINTOUT PRETTY TABLE <<<\n'
470     #header
471     printout = ''
472     printout += '%-25s'%'Process'
473     printout += ':'
474     for item, val in final_histos['nominal'].items():
475 peller 1.70 printout += '%-12s'%item
476 peller 1.69 print printout+'\n'
477 nmohr 1.68 for key in final_histos:
478     printout = ''
479 peller 1.69 printout += '%-25s'%key
480     printout += ':'
481 nmohr 1.68 for item, val in final_histos[key].items():
482 peller 1.74 printout += '%-12s'%str('%0.5f'%val.Integral())
483 nmohr 1.68 print printout
484 peller 1.1
485 nmohr 1.68 #-----------------------------------------------------------------------------------------------------------
486 peller 1.36
487 peller 1.62 # -------------------- write DATAcard: ----------------------------------------------------------------------
488     DCprocessseparatordict = {'WS':':','TH':'/'}
489 nmohr 1.68 # create two datacards: for TH an WS
490 peller 1.62 for DCtype in ['WS','TH']:
491     columns=len(setup)
492     f = open(outpath+'vhbb_DC_%s_%s.txt'%(DCtype,ROOToutname),'w')
493     f.write('imax\t1\tnumber of channels\n')
494     f.write('jmax\t%s\tnumber of backgrounds (\'*\' = automatic)\n'%(columns-1))
495     f.write('kmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
496     f.write('shapes * * vhbb_%s_%s.root $CHANNEL%s$PROCESS $CHANNEL%s$PROCESS$SYSTEMATIC\n\n'%(DCtype,ROOToutname,DCprocessseparatordict[DCtype],DCprocessseparatordict[DCtype]))
497     f.write('bin\t%s\n\n'%Datacardbin)
498 nmohr 1.93 if toy or signal_inject:
499 peller 1.76 f.write('observation\t%s\n\n'%(hDummy.Integral()))
500 peller 1.75 else:
501 peller 1.76 f.write('observation\t%s\n\n'%(theData.Integral()))
502 nmohr 1.68 # datacard bin
503 peller 1.62 f.write('bin')
504     for c in range(0,columns): f.write('\t%s'%Datacardbin)
505     f.write('\n')
506 nmohr 1.68 # datacard process
507 peller 1.62 f.write('process')
508     for c in setup: f.write('\t%s'%Dict[c])
509     f.write('\n')
510     f.write('process')
511     for c in range(0,columns): f.write('\t%s'%c)
512     f.write('\n')
513 nmohr 1.68 # datacard yields
514 peller 1.62 f.write('rate')
515 nmohr 1.68 for c in setup:
516     f.write('\t%s'%final_histos['nominal'][c].Integral())
517 peller 1.29 f.write('\n')
518 nmohr 1.68 # get list of systematics in use
519 peller 1.62 InUse=eval(config.get('Datacard','InUse_%s'%pt_region))
520 nmohr 1.68 # write non-shape systematics
521 peller 1.62 for item in InUse:
522     f.write(item)
523     what=eval(config.get('Datacard',item))
524     f.write('\t%s'%what['type'])
525     for c in setup:
526     if c in what:
527 nmohr 1.92 if '_eff_e' in item and 'Zmm' in data_sample_names: f.write('\t-')
528     elif '_eff_m' in item and 'Zee' in data_sample_names: f.write('\t-')
529     elif '_trigger_e' in item and 'Zmm' in data_sample_names: f.write('\t-')
530     elif '_trigger_m' in item and 'Zee' in data_sample_names: f.write('\t-')
531 peller 1.35 else:
532 peller 1.62 f.write('\t%s'%what[c])
533 peller 1.35 else:
534     f.write('\t-')
535     f.write('\n')
536 peller 1.77 if not ignore_stats:
537 nmohr 1.68 # Write statistical shape variations
538 peller 1.77 if binstat:
539     for c in setup:
540     for bin in range(0,nBins):
541 nmohr 1.91 if bin in binsBelowThreshold[c]:
542     f.write('%s_bin%s_%s_%s\tshape'%(systematicsnaming['stats'],bin,Dict[c],Datacardbin))
543     for it in range(0,columns):
544     if it == setup.index(c):
545     f.write('\t1.0')
546     else:
547     f.write('\t-')
548     f.write('\n')
549 peller 1.77 else:
550     for c in setup:
551 peller 1.87 f.write('%s_%s_%s\tshape'%(systematicsnaming['stats'],Dict[c],Datacardbin))
552 peller 1.62 for it in range(0,columns):
553     if it == setup.index(c):
554     f.write('\t1.0')
555     else:
556     f.write('\t-')
557     f.write('\n')
558 nmohr 1.68 # UEPS systematics
559 nmohr 1.84 for weightF_sys in weightF_systematics:
560     f.write('%s\tshape' %(systematicsnaming[weightF_sys]))
561 peller 1.62 for it in range(0,columns): f.write('\t1.0')
562     f.write('\n')
563 nmohr 1.68 # additional sample systematics
564 peller 1.62 if addSample_sys:
565     alreadyAdded = []
566     for newSample in addSample_sys.iterkeys():
567     for c in setup:
568 peller 1.75 if not c == GroupDict[newSample]: continue
569 peller 1.62 if Dict[c] in alreadyAdded: continue
570 nmohr 1.71 f.write('%s_%s\tshape'%(systematicsnaming['model'],Dict[c]))
571 peller 1.62 for it in range(0,columns):
572     if it == setup.index(c):
573     f.write('\t1.0')
574     else:
575     f.write('\t-')
576     f.write('\n')
577     alreadyAdded.append(Dict[c])
578 nmohr 1.68 # regular systematics
579 peller 1.62 for sys in systematics:
580     sys_factor=sys_factor_dict[sys]
581     f.write('%s\tshape'%systematicsnaming[sys])
582 peller 1.63 for c in setup:
583     if c in sys_affecting[sys]:
584     f.write('\t%s'%sys_factor)
585     else:
586     f.write('\t-')
587 peller 1.62 f.write('\n')
588     f.close()
589 nmohr 1.68 # --------------------------------------------------------------------------
590    
591    
592    
593    
594 nmohr 1.31 outfile.Close()