ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/write_regression_systematics.py
Revision: 1.31
Committed: Wed Jan 16 09:56:56 2013 UTC (12 years, 3 months ago) by peller
Content type: text/x-python
Branch: MAIN
Changes since 1.30: +10 -0 lines
Log Message:
removing obsolete files step*

File Contents

# User Rev Content
1 nmohr 1.1 #!/usr/bin/env python
2     from samplesclass import sample
3     from printcolor import printc
4     import pickle
5     import sys
6     import os
7     import ROOT
8     import math
9     import shutil
10     from array import array
11     import warnings
12 peller 1.15 from optparse import OptionParser
13     from BetterConfigParser import BetterConfigParser
14 nmohr 1.1 warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
15    
16    
17     #usage: ./write_regression_systematic.py path
18    
19    
20     #os.mkdir(path+'/sys')
21 nmohr 1.16 argv = sys.argv
22 peller 1.15 parser = OptionParser()
23 nmohr 1.16 parser.add_option("-P", "--path", dest="path", default="",
24     help="path to samples")
25     parser.add_option("-S", "--samples", dest="names", default="",
26     help="samples you want to run on")
27 nmohr 1.19 parser.add_option("-C", "--config", dest="config", default=[], action="append",
28 peller 1.15 help="configuration defining the plots to make")
29     (opts, args) = parser.parse_args(argv)
30 nmohr 1.16 if opts.config =="":
31 peller 1.15 opts.config = "config"
32     print opts.config
33     config = BetterConfigParser()
34     config.read(opts.config)
35     anaTag = config.get("Analysis","tag")
36 peller 1.23 TrainFlag = eval(config.get('Analysis','TrainFlag'))
37 nmohr 1.16 btagLibrary = config.get('BTagReshaping','library')
38 nmohr 1.28 samplesinfo=config.get('Directories','samplesinfo')
39 nmohr 1.30
40     VHbbNameSpace=config.get('VHbbNameSpace','library')
41     ROOT.gSystem.Load(VHbbNameSpace)
42     AngLikeBkgs=eval(config.get('AngularLike','backgrounds'))
43     ang_yield=eval(config.get('AngularLike','yields'))
44    
45 nmohr 1.16 path=opts.path
46 peller 1.31
47    
48     storagesamples = config.get('Directories','storagesamples')
49    
50    
51 nmohr 1.16 namelist=opts.names.split(',')
52     #load info
53 nmohr 1.28 infofile = open(samplesinfo,'r')
54 nmohr 1.16 info = pickle.load(infofile)
55     infofile.close()
56 nmohr 1.1
57     def deltaPhi(phi1, phi2):
58     result = phi1 - phi2
59     while (result > math.pi): result -= 2*math.pi
60     while (result <= -math.pi): result += 2*math.pi
61     return result
62    
63 peller 1.15 def resolutionBias(eta):
64     if(eta< 1.1): return 0.05
65     if(eta< 2.5): return 0.10
66     if(eta< 5): return 0.30
67     return 0
68    
69     def corrPt(pt,eta,mcPt):
70     return (pt+resolutionBias(math.fabs(eta))*(pt-mcPt))/pt
71    
72 nmohr 1.1 def corrCSV(btag, csv, flav):
73     if(csv < 0.): return csv
74     if(csv > 1.): return csv;
75     if(flav == 0): return csv;
76     if(math.fabs(flav) == 5): return btag.ib.Eval(csv)
77     if(math.fabs(flav) == 4): return btag.ic.Eval(csv)
78     if(math.fabs(flav) != 4 and math.fabs(flav) != 5): return btag.il.Eval(csv)
79     return -10000
80    
81    
82 peller 1.15 def csvReshape(sh, pt, eta, csv, flav):
83     return sh.reshape(float(eta), float(pt), float(csv), int(flav))
84    
85    
86 nmohr 1.1 for job in info:
87 peller 1.15 if not job.name in namelist: continue
88 nmohr 1.1 #print job.name
89 peller 1.15 #if job.name != 'ZH120': continue
90 nmohr 1.1 ROOT.gROOT.ProcessLine(
91     "struct H {\
92     int HiggsFlag;\
93     float mass;\
94     float pt;\
95     float eta;\
96     float phi;\
97     float dR;\
98     float dPhi;\
99     float dEta;\
100     } ;"
101     )
102 peller 1.15 if anaTag == '7TeV':
103 nmohr 1.27 ROOT.gSystem.Load(btagLibrary)
104     from ROOT import BTagShape
105     btagNom = BTagShape("../data/csvdiscr.root")
106     btagNom.computeFunctions()
107     btagUp = BTagShape("../data/csvdiscr.root")
108     btagUp.computeFunctions(+1.,0.)
109     btagDown = BTagShape("../data/csvdiscr.root")
110     btagDown.computeFunctions(-1.,0.)
111     btagFUp = BTagShape("../data/csvdiscr.root")
112     btagFUp.computeFunctions(0.,+1.)
113     btagFDown = BTagShape("../data/csvdiscr.root")
114     btagFDown.computeFunctions(0.,-1.)
115 peller 1.15 elif anaTag == '8TeV':
116 nmohr 1.27 ROOT.gSystem.Load(btagLibrary)
117     from ROOT import BTagShapeInterface
118     btagNom = BTagShapeInterface("../data/csvdiscr.root",0,0)
119     btagUp = BTagShapeInterface("../data/csvdiscr.root",+1,0)
120     btagDown = BTagShapeInterface("../data/csvdiscr.root",-1,0)
121     btagFUp = BTagShapeInterface("../data/csvdiscr.root",0,+1.)
122     btagFDown = BTagShapeInterface("../data/csvdiscr.root",0,-1.)
123 nmohr 1.1
124     print '\t - %s' %(job.name)
125 peller 1.31 <<<<<<< write_regression_systematics.py
126     input = TFile.Open(storagesamples+'/env/'+job.getpath(),'read')
127     output = TFile.Open(path+'/sys/'+job.prefix+job.identifier+'.root','recreate')
128     =======
129 nmohr 1.29 input = ROOT.TFile.Open(path+'/'+job.getpath(),'read')
130     output = ROOT.TFile.Open(path+'/sys/'+job.prefix+job.identifier+'.root','recreate')
131 peller 1.31 >>>>>>> 1.29
132 nmohr 1.1
133     input.cd()
134     obj = ROOT.TObject
135     for key in ROOT.gDirectory.GetListOfKeys():
136     input.cd()
137     obj = key.ReadObj()
138     #print obj.GetName()
139     if obj.GetName() == job.tree:
140     continue
141     output.cd()
142     #print key.GetName()
143     obj.Write(key.GetName())
144    
145 nmohr 1.3 input.cd()
146 nmohr 1.1 tree = input.Get(job.tree)
147     nEntries = tree.GetEntries()
148    
149     job.addpath('/sys')
150     if job.type != 'DATA':
151     job.SYS = ['Nominal','JER_up','JER_down','JES_up','JES_down','beff_up','beff_down','bmis_up','bmis_down']
152    
153     H = ROOT.H()
154     HNoReg = ROOT.H()
155     tree.SetBranchStatus('H',0)
156     output.cd()
157     newtree = tree.CloneTree(0)
158    
159     hJ0 = ROOT.TLorentzVector()
160     hJ1 = ROOT.TLorentzVector()
161    
162 peller 1.15 regWeight = config.get("Regression","regWeight")
163     regDict = eval(config.get("Regression","regDict"))
164     regVars = eval(config.get("Regression","regVars"))
165     useMET = eval(config.get("Regression","useMET"))
166     usePtRaw = eval(config.get("Regression","usePtRaw"))
167 nmohr 1.21 useMt = eval(config.get("Regression","useMt"))
168 peller 1.15 useRho25 = eval(config.get("Regression","useRho25"))
169 nmohr 1.1
170     #Regression branches
171     applyRegression = True
172     hJet_pt = array('f',[0]*2)
173     hJet_e = array('f',[0]*2)
174     newtree.Branch( 'H', H , 'HiggsFlag/I:mass/F:pt/F:eta:phi/F:dR/F:dPhi/F:dEta/F' )
175     newtree.Branch( 'HNoReg', HNoReg , 'HiggsFlag/I:mass/F:pt/F:eta:phi/F:dR/F:dPhi/F:dEta/F' )
176     Event = array('f',[0])
177     METet = array('f',[0])
178     rho25 = array('f',[0])
179     METphi = array('f',[0])
180     fRho25 = ROOT.TTreeFormula("rho25",'rho25',tree)
181     fEvent = ROOT.TTreeFormula("Event",'EVENT.event',tree)
182     fMETet = ROOT.TTreeFormula("METet",'METnoPU.et',tree)
183     fMETphi = ROOT.TTreeFormula("METphi",'METnoPU.phi',tree)
184 nmohr 1.30 fHVMass = ROOT.TTreeFormula("HVMass",'HVMass',tree)
185 nmohr 1.21 hJet_MtArray = [array('f',[0]),array('f',[0])]
186 nmohr 1.24 hJet_EtArray = [array('f',[0]),array('f',[0])]
187 nmohr 1.1 hJet_MET_dPhi = array('f',[0]*2)
188     hJet_regWeight = array('f',[0]*2)
189     hJet_MET_dPhiArray = [array('f',[0]),array('f',[0])]
190 nmohr 1.27 hJet_ptRawArray = [array('f',[0]),array('f',[0])]
191 nmohr 1.1 newtree.Branch('hJet_MET_dPhi',hJet_MET_dPhi,'hJet_MET_dPhi[2]/F')
192     newtree.Branch('hJet_regWeight',hJet_regWeight,'hJet_regWeight[2]/F')
193     readerJet0 = ROOT.TMVA.Reader("!Color:!Silent" )
194     readerJet1 = ROOT.TMVA.Reader("!Color:!Silent" )
195    
196     theForms = {}
197     theVars0 = {}
198     theVars1 = {}
199 nmohr 1.25 def addVarsToReader(reader,theVars,theForms,i,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray):
200 nmohr 1.24 for key in regVars:
201     var = regDict[key]
202     theVars[key] = array( 'f', [ 0 ] )
203     if var == 'Jet_MET_dPhi':
204     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
205     reader.AddVariable( key, hJet_MET_dPhiArray[i] )
206     elif var == 'METet':
207     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
208 nmohr 1.26 reader.AddVariable( key, METet )
209 nmohr 1.24 elif var == 'rho25':
210     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
211 nmohr 1.26 reader.AddVariable( key, rho25 )
212 nmohr 1.24 elif var == 'Jet_mt':
213     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
214 nmohr 1.26 reader.AddVariable( key, hJet_MtArray[i] )
215 nmohr 1.24 elif var == 'Jet_et':
216     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
217 nmohr 1.26 reader.AddVariable( key, hJet_EtArray[i] )
218 nmohr 1.24 elif var == 'Jet_ptRaw':
219     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
220 nmohr 1.26 reader.AddVariable( key, hJet_ptRawArray[i] )
221 nmohr 1.24 else:
222     reader.AddVariable(key,theVars[key])
223     formula = regDict[key].replace("[0]","[%.0f]" %i)
224     print 'Adding var: %s with %s to readerJet%.0f' %(key,formula,i)
225     theForms['form_reg_%s_%.0f'%(key,i)] = ROOT.TTreeFormula("form_reg_%s_%.0f"%(key,i),'%s' %(formula),tree)
226     return
227 nmohr 1.25 addVarsToReader(readerJet0,theVars0,theForms,0,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray)
228     addVarsToReader(readerJet1,theVars1,theForms,1,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray)
229 nmohr 1.24 readerJet0.BookMVA( "jet0Regression", regWeight )
230     readerJet1.BookMVA( "jet1Regression", regWeight )
231    
232    
233 nmohr 1.1
234 bortigno 1.4 #Add training Flag
235 nmohr 1.18 EventForTraining = array('i',[0])
236     newtree.Branch('EventForTraining',EventForTraining,'EventForTraining/I')
237 bortigno 1.4 EventForTraining[0]=0
238 bortigno 1.6
239 bortigno 1.5 TFlag=ROOT.TTreeFormula("EventForTraining","EVENT.event%2",tree)
240 nmohr 1.30
241     #Angular Likelihood
242     angleHB = array('f',[0])
243     newtree.Branch('angleHB',angleHB,'angleHB/F')
244     angleLZ = array('f',[0])
245     newtree.Branch('angleLZ',angleLZ,'angleLZ/F')
246     angleZZS = array('f',[0])
247     newtree.Branch('angleZZS',angleZZS,'angleZZS/F')
248     kinLikeRatio = array('f',[0]*len(AngLikeBkgs))
249     newtree.Branch('kinLikeRatio',kinLikeRatio,'%s/F' %(':'.join(AngLikeBkgs)))
250     fAngleHB = ROOT.TTreeFormula("fAngleHB",'abs(VHbb::ANGLEHB(hJet_pt[0],hJet_eta[0],hJet_phi[0],hJet_e[0]*(hJet_pt[0]/hJet_pt[0]),hJet_pt[1],hJet_eta[1],hJet_phi[1],hJet_e[1]*(hJet_pt[1]/hJet_pt[1])))',newtree)
251     fAngleLZ = ROOT.TTreeFormula("fAngleLZ",'abs(VHbb::ANGLELZ(vLepton_pt[0],vLepton_eta[0],vLepton_phi[0],vLepton_mass[0],vLepton_pt[1],vLepton_eta[1],vLepton_phi[1],vLepton_mass[1]))',newtree)
252     fAngleZZS = ROOT.TTreeFormula("fAngleZZS",'abs(VHbb::ANGLELZ(H.pt,H.eta,H.phi,H.pt,V.pt,V.eta,V.phi,V.mass))',newtree)
253     likeSBH = array('f',[0]*len(AngLikeBkgs))
254     likeBBH = array('f',[0]*len(AngLikeBkgs))
255     likeSLZ = array('f',[0]*len(AngLikeBkgs))
256     likeBLZ = array('f',[0]*len(AngLikeBkgs))
257     likeSZZS = array('f',[0]*len(AngLikeBkgs))
258     likeBZZS = array('f',[0]*len(AngLikeBkgs))
259     likeSMassZS = array('f',[0]*len(AngLikeBkgs))
260     likeBMassZS = array('f',[0]*len(AngLikeBkgs))
261    
262     SigBH = []; BkgBH = []; SigLZ = []; BkgLZ = []; SigZZS = []; BkgZZS = []; SigMassZS = []; BkgMassZS = [];
263     for angLikeBkg in AngLikeBkgs:
264     f = ROOT.TFile("../data/angleFitFunctions-%s.root"%(angLikeBkg),"READ")
265     SigBH.append(f.Get("sigFuncBH"))
266     BkgBH.append(f.Get("bkgFuncBH"))
267     SigLZ.append(f.Get("sigFuncLZ"))
268     BkgLZ.append(f.Get("bkgFuncLZ"))
269     SigZZS.append(f.Get("sigFuncZZS"))
270     BkgZZS.append(f.Get("bkgFuncZZS"))
271     SigMassZS.append(f.Get("sigFuncMassZS"))
272     BkgMassZS.append(f.Get("bkgFuncMassZS"))
273     f.Close()
274    
275 nmohr 1.1
276     if job.type != 'DATA':
277     #CSV branches
278     hJet_flavour = array('f',[0]*2)
279     hJet_csv = array('f',[0]*2)
280     hJet_csvOld = array('f',[0]*2)
281     hJet_csvUp = array('f',[0]*2)
282     hJet_csvDown = array('f',[0]*2)
283     hJet_csvFUp = array('f',[0]*2)
284     hJet_csvFDown = array('f',[0]*2)
285     newtree.Branch('hJet_csvOld',hJet_csvOld,'hJet_csvOld[2]/F')
286     newtree.Branch('hJet_csvUp',hJet_csvUp,'hJet_csvUp[2]/F')
287     newtree.Branch('hJet_csvDown',hJet_csvDown,'hJet_csvDown[2]/F')
288     newtree.Branch('hJet_csvFUp',hJet_csvFUp,'hJet_csvFUp[2]/F')
289     newtree.Branch('hJet_csvFDown',hJet_csvFDown,'hJet_csvFDown[2]/F')
290    
291     #JER branches
292     hJet_pt_JER_up = array('f',[0]*2)
293     newtree.Branch('hJet_pt_JER_up',hJet_pt_JER_up,'hJet_pt_JER_up[2]/F')
294     hJet_pt_JER_down = array('f',[0]*2)
295     newtree.Branch('hJet_pt_JER_down',hJet_pt_JER_down,'hJet_pt_JER_down[2]/F')
296     hJet_e_JER_up = array('f',[0]*2)
297     newtree.Branch('hJet_e_JER_up',hJet_e_JER_up,'hJet_e_JER_up[2]/F')
298     hJet_e_JER_down = array('f',[0]*2)
299     newtree.Branch('hJet_e_JER_down',hJet_e_JER_down,'hJet_e_JER_down[2]/F')
300     H_JER = array('f',[0]*4)
301     newtree.Branch('H_JER',H_JER,'mass_up:mass_down:pt_up:pt_down/F')
302    
303     #JES branches
304     hJet_pt_JES_up = array('f',[0]*2)
305     newtree.Branch('hJet_pt_JES_up',hJet_pt_JES_up,'hJet_pt_JES_up[2]/F')
306     hJet_pt_JES_down = array('f',[0]*2)
307     newtree.Branch('hJet_pt_JES_down',hJet_pt_JES_down,'hJet_pt_JES_down[2]/F')
308     hJet_e_JES_up = array('f',[0]*2)
309     newtree.Branch('hJet_e_JES_up',hJet_e_JES_up,'hJet_e_JES_up[2]/F')
310     hJet_e_JES_down = array('f',[0]*2)
311     newtree.Branch('hJet_e_JES_down',hJet_e_JES_down,'hJet_e_JES_down[2]/F')
312     H_JES = array('f',[0]*4)
313     newtree.Branch('H_JES',H_JES,'mass_up:mass_down:pt_up:pt_down/F')
314 nmohr 1.27 lheWeight = array('f',[0])
315 nmohr 1.17 if job.type != 'DY':
316 nmohr 1.26 newtree.Branch('lheWeight',lheWeight,'lheWeight/F')
317     lheWeight[0] = 1.
318 nmohr 1.1
319    
320     #iter=0
321    
322    
323     for entry in range(0,nEntries):
324     tree.GetEntry(entry)
325    
326 peller 1.22 if job.type != 'DATA' and TrainFlag:
327 nmohr 1.18 EventForTraining[0]=int(not TFlag.EvalInstance())
328 bortigno 1.9
329 nmohr 1.1
330     #get
331     hJet_pt = tree.hJet_pt
332     hJet_e = tree.hJet_e
333 nmohr 1.27 hJet_pt1 = tree.hJet_pt[1]
334 nmohr 1.1 hJet_pt0 = tree.hJet_pt[0]
335     hJet_pt1 = tree.hJet_pt[1]
336     hJet_eta0 = tree.hJet_eta[0]
337     hJet_eta1 = tree.hJet_eta[1]
338     hJet_genPt0 = tree.hJet_genPt[0]
339     hJet_genPt1 = tree.hJet_genPt[1]
340 nmohr 1.27 hJet_ptRaw0 = tree.hJet_ptRaw[0]
341     hJet_ptRaw1 = tree.hJet_ptRaw[1]
342 nmohr 1.1 hJet_e0 = tree.hJet_e[0]
343     hJet_e1 = tree.hJet_e[1]
344     hJet_phi0 = tree.hJet_phi[0]
345     hJet_phi1 = tree.hJet_phi[1]
346     hJet_JECUnc0 = tree.hJet_JECUnc[0]
347     hJet_JECUnc1 = tree.hJet_JECUnc[1]
348    
349     Event[0]=fEvent.EvalInstance()
350     METet[0]=fMETet.EvalInstance()
351     rho25[0]=fRho25.EvalInstance()
352     METphi[0]=fMETphi.EvalInstance()
353     for key, value in regDict.items():
354 nmohr 1.26 if not (value == 'Jet_MET_dPhi' or value == 'METet' or value == "rho25" or value == "Jet_et" or value == 'Jet_mt' or value == 'Jet_ptRaw'):
355 nmohr 1.24 theVars0[key][0] = theForms["form_reg_%s_0" %(key)].EvalInstance()
356     theVars1[key][0] = theForms["form_reg_%s_1" %(key)].EvalInstance()
357 nmohr 1.27 hJet_MET_dPhi[0] = deltaPhi(METphi[0],hJet_phi0)
358     hJet_MET_dPhi[1] = deltaPhi(METphi[0],hJet_phi1)
359     hJet_MET_dPhiArray[0][0] = deltaPhi(METphi[0],hJet_phi0)
360     hJet_MET_dPhiArray[1][0] = deltaPhi(METphi[0],hJet_phi1)
361     if not job.type == 'DATA':
362     corrRes0 = corrPt(hJet_pt0,hJet_eta0,hJet_genPt0)
363     corrRes1 = corrPt(hJet_pt1,hJet_eta1,hJet_genPt1)
364     hJet_ptRaw0 *= corrRes0
365     hJet_ptRaw1 *= corrRes1
366     hJet_ptRawArray[0][0] = hJet_ptRaw0
367     hJet_ptRawArray[1][0] = hJet_ptRaw1
368 nmohr 1.1
369 peller 1.15
370 nmohr 1.1 if applyRegression:
371     hJ0.SetPtEtaPhiE(hJet_pt0,hJet_eta0,hJet_phi0,hJet_e0)
372     hJ1.SetPtEtaPhiE(hJet_pt1,hJet_eta1,hJet_phi1,hJet_e1)
373     HNoReg.HiggsFlag = 1
374     HNoReg.mass = (hJ0+hJ1).M()
375     HNoReg.pt = (hJ0+hJ1).Pt()
376     HNoReg.eta = (hJ0+hJ1).Eta()
377     HNoReg.phi = (hJ0+hJ1).Phi()
378     HNoReg.dR = hJ0.DeltaR(hJ1)
379     HNoReg.dPhi = hJ0.DeltaPhi(hJ1)
380     HNoReg.dEta = abs(hJ0.Eta()-hJ1.Eta())
381 nmohr 1.21 hJet_MtArray[0][0] = hJ0.Mt()
382     hJet_MtArray[1][0] = hJ1.Mt()
383 nmohr 1.24 hJet_EtArray[0][0] = hJ0.Et()
384     hJet_EtArray[1][0] = hJ1.Et()
385 peller 1.15 rPt0 = max(0.0001,readerJet0.EvaluateRegression( "jet0Regression" )[0])
386     rPt1 = max(0.0001,readerJet1.EvaluateRegression( "jet1Regression" )[0])
387 nmohr 1.1 hJet_regWeight[0] = rPt0/hJet_pt0
388     hJet_regWeight[1] = rPt1/hJet_pt1
389     rE0 = hJet_e0*hJet_regWeight[0]
390     rE1 = hJet_e1*hJet_regWeight[1]
391     hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
392     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
393     tree.hJet_pt[0] = rPt0
394     tree.hJet_pt[1] = rPt1
395     tree.hJet_e[0] = rE0
396     tree.hJet_e[1] = rE1
397     H.HiggsFlag = 1
398     H.mass = (hJ0+hJ1).M()
399     H.pt = (hJ0+hJ1).Pt()
400     H.eta = (hJ0+hJ1).Eta()
401     H.phi = (hJ0+hJ1).Phi()
402     H.dR = hJ0.DeltaR(hJ1)
403     H.dPhi = hJ0.DeltaPhi(hJ1)
404     H.dEta = abs(hJ0.Eta()-hJ1.Eta())
405     if hJet_regWeight[0] > 5. or hJet_regWeight[1] > 5.:
406 nmohr 1.24 print 'Event %.0f' %(Event[0])
407 nmohr 1.1 print 'MET %.2f' %(METet[0])
408     print 'rho25 %.2f' %(rho25[0])
409     for key, value in regDict.items():
410 nmohr 1.26 if not (value == 'Jet_MET_dPhi' or value == 'METet' or value == "rho25" or value == "Jet_et" or value == 'Jet_mt' or value == 'Jet_ptRaw'):
411 nmohr 1.24 print '%s 0: %.2f'%(key, theVars0[key][0])
412 nmohr 1.26 print '%s 1: %.2f'%(key, theVars1[key][0])
413 nmohr 1.1 for i in range(2):
414     print 'dPhi %.0f %.2f' %(i,hJet_MET_dPhiArray[i][0])
415 nmohr 1.24 for i in range(2):
416     print 'ptRaw %.0f %.2f' %(i,hJet_ptRawArray[i][0])
417     for i in range(2):
418     print 'Mt %.0f %.2f' %(i,hJet_MtArray[i][0])
419     for i in range(2):
420     print 'Et %.0f %.2f' %(i,hJet_EtArray[i][0])
421 nmohr 1.1 print 'corr 0 %.2f' %(hJet_regWeight[0])
422     print 'corr 1 %.2f' %(hJet_regWeight[1])
423     print 'rPt0 %.2f' %(rPt0)
424     print 'rPt1 %.2f' %(rPt1)
425     print 'rE0 %.2f' %(rE0)
426     print 'rE1 %.2f' %(rE1)
427     print 'Mass %.2f' %(H.mass)
428 nmohr 1.30
429     angleHB[0]=fAngleHB.EvalInstance()
430     angleLZ[0]=fAngleLZ.EvalInstance()
431     angleZZS[0]=fAngleZZS.EvalInstance()
432    
433     for i, angLikeBkg in enumerate(AngLikeBkgs):
434     likeSBH[i] = math.fabs(SigBH[i].Eval(angleHB[0]))
435     likeBBH[i] = math.fabs(BkgBH[i].Eval(angleHB[0]))
436    
437     likeSZZS[i] = math.fabs(SigZZS[i].Eval(angleZZS[0]))
438     likeBZZS[i] = math.fabs(BkgZZS[i].Eval(angleZZS[0]))
439    
440     likeSLZ[i] = math.fabs(SigLZ[i].Eval(angleLZ[0]))
441     likeBLZ[i] = math.fabs(BkgLZ[i].Eval(angleLZ[0]))
442    
443     likeSMassZS[i] = math.fabs(SigMassZS[i].Eval(fHVMass.EvalInstance()))
444     likeBMassZS[i] = math.fabs(BkgMassZS[i].Eval(fHVMass.EvalInstance()))
445    
446     scaleSig = float( ang_yield['Sig'] / (ang_yield['Sig'] + ang_yield[angLikeBkg]))
447     scaleBkg = float( ang_yield[angLikeBkg] / (ang_yield['Sig'] + ang_yield[angLikeBkg]) )
448    
449     numerator = (likeSBH[i]*likeSZZS[i]*likeSLZ[i]*likeSMassZS[i]);
450     denominator = ((scaleBkg*likeBLZ[i]*likeBZZS[i]*likeBBH[i]*likeBMassZS[i])+(scaleSig*likeSBH[i]*likeSZZS[i]*likeSLZ[i]*likeSMassZS[i]))
451    
452     if denominator > 0:
453     kinLikeRatio[i] = numerator/denominator;
454     else:
455     kinLikeRatio[i] = 0;
456    
457 nmohr 1.1 if job.type == 'DATA':
458     newtree.Fill()
459     continue
460    
461     for i in range(2):
462 peller 1.15 flavour = int(tree.hJet_flavour[i])
463     pt = float(tree.hJet_pt[i])
464     eta = float(tree.hJet_eta[i])
465     csv = float(tree.hJet_csv[i])
466 nmohr 1.26 hJet_csvOld[i] = csv
467     if anaTag == '7TeV':
468     tree.hJet_csv[i] = corrCSV(btagNom,csv,flavour)
469     hJet_csvDown[i] = corrCSV(btagDown,csv,flavour)
470     hJet_csvUp[i] = corrCSV(btagUp,csv,flavour)
471     hJet_csvFDown[i] = corrCSV(btagFDown,csv,flavour)
472     hJet_csvFUp[i] = corrCSV(btagFUp,csv,flavour)
473     elif anaTag == '8TeV':
474     tree.hJet_csv[i] = btagNom.reshape(eta,pt,csv,flavour)
475     hJet_csvDown[i] = btagDown.reshape(eta,pt,csv,flavour)
476     hJet_csvUp[i] = btagUp.reshape(eta,pt,csv,flavour)
477     hJet_csvFDown[i] = btagFDown.reshape(eta,pt,csv,flavour)
478     hJet_csvFUp[i] = btagFUp.reshape(eta,pt,csv,flavour)
479 nmohr 1.1
480     for updown in ['up','down']:
481     #JER
482     if updown == 'up':
483     inner = 0.06
484     outer = 0.1
485     if updown == 'down':
486     inner = -0.06
487     outer = -0.1
488     #Calculate
489     if abs(hJet_eta0)<1.1: res0 = inner
490     else: res0 = outer
491     if abs(hJet_eta1)<1.1: res1 = inner
492     else: res1 = outer
493     rPt0 = hJet_pt0 + (hJet_pt0-hJet_genPt0)*res0
494     rPt1 = hJet_pt1 + (hJet_pt1-hJet_genPt1)*res1
495     rE0 = hJet_e0*rPt0/hJet_pt0
496     rE1 = hJet_e1*rPt1/hJet_pt1
497 nmohr 1.27 if applyRegression:
498     for key in regVars:
499     var = regDict[key]
500     if var == 'Jet_pt' or var == 'Jet_e' or var == 'hJet_pt' or var == 'hJet_e' or var == 'Jet_ptRaw':
501     if var == 'Jet_ptRaw':
502     hJet_ptRawArray[0][0] = hJet_ptRaw0*corrRes0*rPt0/hJet_pt0
503     hJet_ptRawArray[1][0] = hJet_ptRaw1*rPt1/hJet_pt1
504    
505     elif var == 'Jet_pt' or var == 'hJet_pt':
506     theVars0[key] = rPt0
507     theVars1[key] = rPt1
508     elif var == 'Jet_e' or var == 'hJet_e':
509     theVars0[key] = rE0
510     theVars1[key] = rE1
511     hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
512     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
513     hJet_MtArray[0][0] = hJ0.Mt()
514     hJet_MtArray[1][0] = hJ1.Mt()
515     hJet_EtArray[0][0] = hJ0.Et()
516     hJet_EtArray[1][0] = hJ1.Et()
517     rPt0 = max(0.0001,readerJet0.EvaluateRegression( "jet0Regression" )[0])
518     rPt1 = max(0.0001,readerJet1.EvaluateRegression( "jet1Regression" )[0])
519     rE0 = hJet_e0*rPt0/hJet_pt0
520     rE1 = hJet_e1*rPt1/hJet_pt1
521 nmohr 1.1 hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
522     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
523     #Set
524     if updown == 'up':
525     hJet_pt_JER_up[0]=rPt0
526     hJet_pt_JER_up[1]=rPt1
527     hJet_e_JER_up[0]=rE0
528     hJet_e_JER_up[1]=rE1
529     H_JER[0]=(hJ0+hJ1).M()
530     H_JER[2]=(hJ0+hJ1).Pt()
531     if updown == 'down':
532     hJet_pt_JER_down[0]=rPt0
533     hJet_pt_JER_down[1]=rPt1
534     hJet_e_JER_down[0]=rE0
535     hJet_e_JER_down[1]=rE1
536     H_JER[1]=(hJ0+hJ1).M()
537     H_JER[3]=(hJ0+hJ1).Pt()
538    
539     #JES
540     if updown == 'up':
541     variation=1
542     if updown == 'down':
543     variation=-1
544     #calculate
545     rPt0 = hJet_pt0*(1+variation*hJet_JECUnc0)
546     rPt1 = hJet_pt1*(1+variation*hJet_JECUnc1)
547     rE0 = hJet_e0*(1+variation*hJet_JECUnc0)
548     rE1 = hJet_e1*(1+variation*hJet_JECUnc1)
549 nmohr 1.27 if applyRegression:
550     for key in regVars:
551     var = regDict[key]
552     if var == 'Jet_pt' or var == 'Jet_e' or var == 'hJet_pt' or var == 'hJet_e' or var == 'Jet_ptRaw':
553     if var == 'Jet_ptRaw':
554     hJet_ptRawArray[0][0] = hJet_ptRaw0*(1+variation*hJet_JECUnc0)
555     hJet_ptRawArray[1][0] = hJet_ptRaw1*(1+variation*hJet_JECUnc1)
556    
557     elif var == 'Jet_pt' or var == 'hJet_pt':
558     theVars0[key] = rPt0
559     theVars1[key] = rPt1
560     elif var == 'Jet_e' or var == 'hJet_e':
561     theVars0[key] = rE0
562     theVars1[key] = rE1
563     hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
564     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
565     hJet_MtArray[0][0] = hJ0.Mt()
566     hJet_MtArray[1][0] = hJ1.Mt()
567     hJet_EtArray[0][0] = hJ0.Et()
568     hJet_EtArray[1][0] = hJ1.Et()
569     rPt0 = max(0.0001,readerJet0.EvaluateRegression( "jet0Regression" )[0])
570     rPt1 = max(0.0001,readerJet1.EvaluateRegression( "jet1Regression" )[0])
571     rE0 = hJet_e0*rPt0/hJet_pt0
572     rE1 = hJet_e1*rPt1/hJet_pt1
573 nmohr 1.1 hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
574     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
575     #Fill
576     if updown == 'up':
577     hJet_pt_JES_up[0]=rPt0
578     hJet_pt_JES_up[1]=rPt1
579     hJet_e_JES_up[0]=rE0
580     hJet_e_JES_up[1]=rE1
581     H_JES[0]=(hJ0+hJ1).M()
582     H_JES[2]=(hJ0+hJ1).Pt()
583     if updown == 'down':
584     hJet_pt_JES_down[0]=rPt0
585     hJet_pt_JES_down[1]=rPt1
586     hJet_e_JES_down[0]=rE0
587     hJet_e_JES_down[1]=rE1
588     H_JES[1]=(hJ0+hJ1).M()
589     H_JES[3]=(hJ0+hJ1).Pt()
590    
591     newtree.Fill()
592    
593     newtree.AutoSave()
594     output.Close()
595    
596     #dump info
597 nmohr 1.28 #infofile = open(samplesinfo,'w')
598 peller 1.15 #pickle.dump(info,infofile)
599     #infofile.close()