ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/write_regression_systematics.py
Revision: 1.36
Committed: Mon Mar 4 14:35:44 2013 UTC (12 years, 2 months ago) by nmohr
Content type: text/x-python
Branch: MAIN
Changes since 1.35: +89 -65 lines
Log Message:
Add lheweight in sys, to be validated

File Contents

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