ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/write_regression_systematics.py
Revision: 1.34
Committed: Fri Jan 25 16:18:00 2013 UTC (12 years, 3 months ago) by nmohr
Content type: text/x-python
Branch: MAIN
Changes since 1.33: +66 -18 lines
Log Message:
Restructuring, still to be validated, workspace writing missing

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 peller 1.15 from optparse import OptionParser
11 nmohr 1.34 ROOT.gROOT.SetBatch(True)
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 #print job.name
91 peller 1.15 #if job.name != 'ZH120': continue
92 nmohr 1.1 ROOT.gROOT.ProcessLine(
93     "struct H {\
94     int HiggsFlag;\
95     float mass;\
96     float pt;\
97     float eta;\
98     float phi;\
99     float dR;\
100     float dPhi;\
101     float dEta;\
102     } ;"
103     )
104 peller 1.15 if anaTag == '7TeV':
105 nmohr 1.27 ROOT.gSystem.Load(btagLibrary)
106     from ROOT import BTagShape
107     btagNom = BTagShape("../data/csvdiscr.root")
108     btagNom.computeFunctions()
109     btagUp = BTagShape("../data/csvdiscr.root")
110     btagUp.computeFunctions(+1.,0.)
111     btagDown = BTagShape("../data/csvdiscr.root")
112     btagDown.computeFunctions(-1.,0.)
113     btagFUp = BTagShape("../data/csvdiscr.root")
114     btagFUp.computeFunctions(0.,+1.)
115     btagFDown = BTagShape("../data/csvdiscr.root")
116     btagFDown.computeFunctions(0.,-1.)
117 peller 1.15 elif anaTag == '8TeV':
118 nmohr 1.27 ROOT.gSystem.Load(btagLibrary)
119     from ROOT import BTagShapeInterface
120     btagNom = BTagShapeInterface("../data/csvdiscr.root",0,0)
121     btagUp = BTagShapeInterface("../data/csvdiscr.root",+1,0)
122     btagDown = BTagShapeInterface("../data/csvdiscr.root",-1,0)
123     btagFUp = BTagShapeInterface("../data/csvdiscr.root",0,+1.)
124     btagFDown = BTagShapeInterface("../data/csvdiscr.root",0,-1.)
125 nmohr 1.1
126     print '\t - %s' %(job.name)
127 nmohr 1.34 input = ROOT.TFile.Open(pathIN+job.get_path,'read')
128     output = ROOT.TFile.Open(pathOUT+job.get_path+'.root','recreate')
129 peller 1.33 #input = ROOT.TFile.Open(storagesamples+'/env/'+job.getpath(),'read')
130     #output = ROOT.TFile.Open(path+'/sys/'+job.prefix+job.identifier+'.root','recreate')
131 nmohr 1.1
132     input.cd()
133     obj = ROOT.TObject
134     for key in ROOT.gDirectory.GetListOfKeys():
135     input.cd()
136     obj = key.ReadObj()
137     #print obj.GetName()
138     if obj.GetName() == job.tree:
139     continue
140     output.cd()
141     #print key.GetName()
142     obj.Write(key.GetName())
143    
144 nmohr 1.3 input.cd()
145 nmohr 1.1 tree = input.Get(job.tree)
146     nEntries = tree.GetEntries()
147    
148     H = ROOT.H()
149     HNoReg = ROOT.H()
150     tree.SetBranchStatus('H',0)
151     output.cd()
152     newtree = tree.CloneTree(0)
153    
154     hJ0 = ROOT.TLorentzVector()
155     hJ1 = ROOT.TLorentzVector()
156 nmohr 1.34 hFJ0 = ROOT.TLorentzVector()
157     hFJ1 = ROOT.TLorentzVector()
158 nmohr 1.1
159 peller 1.15 regWeight = config.get("Regression","regWeight")
160     regDict = eval(config.get("Regression","regDict"))
161     regVars = eval(config.get("Regression","regVars"))
162 nmohr 1.34 regWeightFilterJets = config.get("Regression","regWeightFilterJets")
163     regDictFilterJets = eval(config.get("Regression","regDictFilterJets"))
164     regVarsFilterJets = eval(config.get("Regression","regVarsFilterJets"))
165 nmohr 1.1
166     #Regression branches
167     applyRegression = True
168     hJet_pt = array('f',[0]*2)
169     hJet_e = array('f',[0]*2)
170 nmohr 1.34 newtree.Branch( 'H', H , 'HiggsFlag/I:mass/F:pt/F:eta/F:phi/F:dR/F:dPhi/F:dEta/F' )
171     newtree.Branch( 'HNoReg', HNoReg , 'HiggsFlag/I:mass/F:pt/F:eta/F:phi/F:dR/F:dPhi/F:dEta/F' )
172     FatHReg = array('f',[0]*2)
173     newtree.Branch('FatHReg',FatHReg,'filteredmass:filteredpt/F')
174 nmohr 1.1 Event = array('f',[0])
175     METet = array('f',[0])
176     rho25 = array('f',[0])
177     METphi = array('f',[0])
178     fRho25 = ROOT.TTreeFormula("rho25",'rho25',tree)
179     fEvent = ROOT.TTreeFormula("Event",'EVENT.event',tree)
180 nmohr 1.34 fFatHFlag = ROOT.TTreeFormula("FatHFlag",'FatH.FatHiggsFlag',tree)
181     fFatHnFilterJets = ROOT.TTreeFormula("FatHnFilterJets",'nfathFilterJets',tree)
182 nmohr 1.1 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 nmohr 1.34 fathFilterJets_regWeight = array('f',[0]*2)
190 nmohr 1.1 hJet_MET_dPhiArray = [array('f',[0]),array('f',[0])]
191 nmohr 1.27 hJet_ptRawArray = [array('f',[0]),array('f',[0])]
192 nmohr 1.1 newtree.Branch('hJet_MET_dPhi',hJet_MET_dPhi,'hJet_MET_dPhi[2]/F')
193     newtree.Branch('hJet_regWeight',hJet_regWeight,'hJet_regWeight[2]/F')
194     readerJet0 = ROOT.TMVA.Reader("!Color:!Silent" )
195     readerJet1 = ROOT.TMVA.Reader("!Color:!Silent" )
196    
197     theForms = {}
198     theVars0 = {}
199     theVars1 = {}
200 nmohr 1.34 def addVarsToReader(reader,regDict,regVars,theVars,theForms,i,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray):
201 nmohr 1.24 for key in regVars:
202     var = regDict[key]
203     theVars[key] = array( 'f', [ 0 ] )
204     if var == 'Jet_MET_dPhi':
205     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
206     reader.AddVariable( key, hJet_MET_dPhiArray[i] )
207     elif var == 'METet':
208     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
209 nmohr 1.26 reader.AddVariable( key, METet )
210 nmohr 1.24 elif var == 'rho25':
211     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
212 nmohr 1.26 reader.AddVariable( key, rho25 )
213 nmohr 1.24 elif var == 'Jet_mt':
214     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
215 nmohr 1.26 reader.AddVariable( key, hJet_MtArray[i] )
216 nmohr 1.24 elif var == 'Jet_et':
217     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
218 nmohr 1.26 reader.AddVariable( key, hJet_EtArray[i] )
219 nmohr 1.24 elif var == 'Jet_ptRaw':
220     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
221 nmohr 1.26 reader.AddVariable( key, hJet_ptRawArray[i] )
222 nmohr 1.24 else:
223     reader.AddVariable(key,theVars[key])
224     formula = regDict[key].replace("[0]","[%.0f]" %i)
225     print 'Adding var: %s with %s to readerJet%.0f' %(key,formula,i)
226     theForms['form_reg_%s_%.0f'%(key,i)] = ROOT.TTreeFormula("form_reg_%s_%.0f"%(key,i),'%s' %(formula),tree)
227     return
228 nmohr 1.34 addVarsToReader(readerJet0,regDict,regVars,theVars0,theForms,0,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray)
229     addVarsToReader(readerJet1,regDict,regVars,theVars1,theForms,1,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray)
230 nmohr 1.24 readerJet0.BookMVA( "jet0Regression", regWeight )
231     readerJet1.BookMVA( "jet1Regression", regWeight )
232    
233 nmohr 1.34 readerFJ0 = ROOT.TMVA.Reader("!Color:!Silent" )
234     readerFJ1 = ROOT.TMVA.Reader("!Color:!Silent" )
235     theFormsFJ = {}
236     theVars0FJ = {}
237     theVars1FJ = {}
238     addVarsToReader(readerFJ0,regDictFilterJets,regVarsFilterJets,theVars0FJ,theFormsFJ,0,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray)
239     addVarsToReader(readerFJ1,regDictFilterJets,regVarsFilterJets,theVars1FJ,theFormsFJ,1,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray)
240     readerFJ0.BookMVA( "jet0RegressionFJ", regWeightFilterJets )
241     readerFJ1.BookMVA( "jet1RegressionFJ", regWeightFilterJets )
242 nmohr 1.24
243 nmohr 1.1
244 bortigno 1.4 #Add training Flag
245 nmohr 1.18 EventForTraining = array('i',[0])
246     newtree.Branch('EventForTraining',EventForTraining,'EventForTraining/I')
247 bortigno 1.4 EventForTraining[0]=0
248 bortigno 1.6
249 bortigno 1.5 TFlag=ROOT.TTreeFormula("EventForTraining","EVENT.event%2",tree)
250 nmohr 1.30
251     #Angular Likelihood
252     angleHB = array('f',[0])
253     newtree.Branch('angleHB',angleHB,'angleHB/F')
254     angleLZ = array('f',[0])
255     newtree.Branch('angleLZ',angleLZ,'angleLZ/F')
256     angleZZS = array('f',[0])
257     newtree.Branch('angleZZS',angleZZS,'angleZZS/F')
258     kinLikeRatio = array('f',[0]*len(AngLikeBkgs))
259     newtree.Branch('kinLikeRatio',kinLikeRatio,'%s/F' %(':'.join(AngLikeBkgs)))
260     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)
261     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)
262     fAngleZZS = ROOT.TTreeFormula("fAngleZZS",'abs(VHbb::ANGLELZ(H.pt,H.eta,H.phi,H.pt,V.pt,V.eta,V.phi,V.mass))',newtree)
263     likeSBH = array('f',[0]*len(AngLikeBkgs))
264     likeBBH = array('f',[0]*len(AngLikeBkgs))
265     likeSLZ = array('f',[0]*len(AngLikeBkgs))
266     likeBLZ = array('f',[0]*len(AngLikeBkgs))
267     likeSZZS = array('f',[0]*len(AngLikeBkgs))
268     likeBZZS = array('f',[0]*len(AngLikeBkgs))
269     likeSMassZS = array('f',[0]*len(AngLikeBkgs))
270     likeBMassZS = array('f',[0]*len(AngLikeBkgs))
271    
272     SigBH = []; BkgBH = []; SigLZ = []; BkgLZ = []; SigZZS = []; BkgZZS = []; SigMassZS = []; BkgMassZS = [];
273     for angLikeBkg in AngLikeBkgs:
274     f = ROOT.TFile("../data/angleFitFunctions-%s.root"%(angLikeBkg),"READ")
275     SigBH.append(f.Get("sigFuncBH"))
276     BkgBH.append(f.Get("bkgFuncBH"))
277     SigLZ.append(f.Get("sigFuncLZ"))
278     BkgLZ.append(f.Get("bkgFuncLZ"))
279     SigZZS.append(f.Get("sigFuncZZS"))
280     BkgZZS.append(f.Get("bkgFuncZZS"))
281     SigMassZS.append(f.Get("sigFuncMassZS"))
282     BkgMassZS.append(f.Get("bkgFuncMassZS"))
283     f.Close()
284    
285 nmohr 1.1
286     if job.type != 'DATA':
287     #CSV branches
288     hJet_flavour = array('f',[0]*2)
289     hJet_csv = array('f',[0]*2)
290     hJet_csvOld = array('f',[0]*2)
291     hJet_csvUp = array('f',[0]*2)
292     hJet_csvDown = array('f',[0]*2)
293     hJet_csvFUp = array('f',[0]*2)
294     hJet_csvFDown = array('f',[0]*2)
295     newtree.Branch('hJet_csvOld',hJet_csvOld,'hJet_csvOld[2]/F')
296     newtree.Branch('hJet_csvUp',hJet_csvUp,'hJet_csvUp[2]/F')
297     newtree.Branch('hJet_csvDown',hJet_csvDown,'hJet_csvDown[2]/F')
298     newtree.Branch('hJet_csvFUp',hJet_csvFUp,'hJet_csvFUp[2]/F')
299     newtree.Branch('hJet_csvFDown',hJet_csvFDown,'hJet_csvFDown[2]/F')
300    
301     #JER branches
302     hJet_pt_JER_up = array('f',[0]*2)
303     newtree.Branch('hJet_pt_JER_up',hJet_pt_JER_up,'hJet_pt_JER_up[2]/F')
304     hJet_pt_JER_down = array('f',[0]*2)
305     newtree.Branch('hJet_pt_JER_down',hJet_pt_JER_down,'hJet_pt_JER_down[2]/F')
306     hJet_e_JER_up = array('f',[0]*2)
307     newtree.Branch('hJet_e_JER_up',hJet_e_JER_up,'hJet_e_JER_up[2]/F')
308     hJet_e_JER_down = array('f',[0]*2)
309     newtree.Branch('hJet_e_JER_down',hJet_e_JER_down,'hJet_e_JER_down[2]/F')
310     H_JER = array('f',[0]*4)
311     newtree.Branch('H_JER',H_JER,'mass_up:mass_down:pt_up:pt_down/F')
312    
313     #JES branches
314     hJet_pt_JES_up = array('f',[0]*2)
315     newtree.Branch('hJet_pt_JES_up',hJet_pt_JES_up,'hJet_pt_JES_up[2]/F')
316     hJet_pt_JES_down = array('f',[0]*2)
317     newtree.Branch('hJet_pt_JES_down',hJet_pt_JES_down,'hJet_pt_JES_down[2]/F')
318     hJet_e_JES_up = array('f',[0]*2)
319     newtree.Branch('hJet_e_JES_up',hJet_e_JES_up,'hJet_e_JES_up[2]/F')
320     hJet_e_JES_down = array('f',[0]*2)
321     newtree.Branch('hJet_e_JES_down',hJet_e_JES_down,'hJet_e_JES_down[2]/F')
322     H_JES = array('f',[0]*4)
323     newtree.Branch('H_JES',H_JES,'mass_up:mass_down:pt_up:pt_down/F')
324 nmohr 1.27 lheWeight = array('f',[0])
325 nmohr 1.17 if job.type != 'DY':
326 nmohr 1.26 newtree.Branch('lheWeight',lheWeight,'lheWeight/F')
327     lheWeight[0] = 1.
328 nmohr 1.1
329    
330     #iter=0
331    
332    
333     for entry in range(0,nEntries):
334     tree.GetEntry(entry)
335    
336 peller 1.22 if job.type != 'DATA' and TrainFlag:
337 nmohr 1.18 EventForTraining[0]=int(not TFlag.EvalInstance())
338 bortigno 1.9
339 nmohr 1.34 #Has fat higgs
340     fatHiggsFlag=fFatHFlag.EvalInstance()*fFatHnFilterJets.EvalInstance()
341 nmohr 1.1
342     #get
343     hJet_pt = tree.hJet_pt
344     hJet_e = tree.hJet_e
345     hJet_pt0 = tree.hJet_pt[0]
346     hJet_pt1 = tree.hJet_pt[1]
347     hJet_eta0 = tree.hJet_eta[0]
348     hJet_eta1 = tree.hJet_eta[1]
349     hJet_genPt0 = tree.hJet_genPt[0]
350     hJet_genPt1 = tree.hJet_genPt[1]
351 nmohr 1.27 hJet_ptRaw0 = tree.hJet_ptRaw[0]
352     hJet_ptRaw1 = tree.hJet_ptRaw[1]
353 nmohr 1.1 hJet_e0 = tree.hJet_e[0]
354     hJet_e1 = tree.hJet_e[1]
355     hJet_phi0 = tree.hJet_phi[0]
356     hJet_phi1 = tree.hJet_phi[1]
357     hJet_JECUnc0 = tree.hJet_JECUnc[0]
358     hJet_JECUnc1 = tree.hJet_JECUnc[1]
359 nmohr 1.34 #Filterjets
360     if fatHiggsFlag:
361     fathFilterJets_pt0 = tree.fathFilterJets_pt[0]
362     fathFilterJets_pt1 = tree.fathFilterJets_pt[1]
363     fathFilterJets_eta0 = tree.fathFilterJets_eta[0]
364     fathFilterJets_eta1 = tree.fathFilterJets_eta[1]
365     fathFilterJets_phi0 = tree.fathFilterJets_phi[0]
366     fathFilterJets_phi1 = tree.fathFilterJets_phi[1]
367     fathFilterJets_e0 = tree.fathFilterJets_e[0]
368     fathFilterJets_e1 = tree.fathFilterJets_e[1]
369 nmohr 1.1
370     Event[0]=fEvent.EvalInstance()
371     METet[0]=fMETet.EvalInstance()
372     rho25[0]=fRho25.EvalInstance()
373     METphi[0]=fMETphi.EvalInstance()
374     for key, value in regDict.items():
375 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'):
376 nmohr 1.24 theVars0[key][0] = theForms["form_reg_%s_0" %(key)].EvalInstance()
377     theVars1[key][0] = theForms["form_reg_%s_1" %(key)].EvalInstance()
378 nmohr 1.34 for key, value in regDictFilterJets.items():
379     if not (value == 'Jet_MET_dPhi' or value == 'METet' or value == "rho25" or value == "Jet_et" or value == 'Jet_mt' or value == 'Jet_ptRaw'):
380     theVars0FJ[key][0] = theFormsFJ["form_reg_%s_0" %(key)].EvalInstance()
381     theVars1FJ[key][0] = theFormsFJ["form_reg_%s_1" %(key)].EvalInstance()
382 nmohr 1.27 hJet_MET_dPhi[0] = deltaPhi(METphi[0],hJet_phi0)
383     hJet_MET_dPhi[1] = deltaPhi(METphi[0],hJet_phi1)
384     hJet_MET_dPhiArray[0][0] = deltaPhi(METphi[0],hJet_phi0)
385     hJet_MET_dPhiArray[1][0] = deltaPhi(METphi[0],hJet_phi1)
386     if not job.type == 'DATA':
387     corrRes0 = corrPt(hJet_pt0,hJet_eta0,hJet_genPt0)
388     corrRes1 = corrPt(hJet_pt1,hJet_eta1,hJet_genPt1)
389     hJet_ptRaw0 *= corrRes0
390     hJet_ptRaw1 *= corrRes1
391     hJet_ptRawArray[0][0] = hJet_ptRaw0
392     hJet_ptRawArray[1][0] = hJet_ptRaw1
393 nmohr 1.1
394 peller 1.15
395 nmohr 1.1 if applyRegression:
396     hJ0.SetPtEtaPhiE(hJet_pt0,hJet_eta0,hJet_phi0,hJet_e0)
397     hJ1.SetPtEtaPhiE(hJet_pt1,hJet_eta1,hJet_phi1,hJet_e1)
398     HNoReg.HiggsFlag = 1
399     HNoReg.mass = (hJ0+hJ1).M()
400     HNoReg.pt = (hJ0+hJ1).Pt()
401     HNoReg.eta = (hJ0+hJ1).Eta()
402     HNoReg.phi = (hJ0+hJ1).Phi()
403     HNoReg.dR = hJ0.DeltaR(hJ1)
404     HNoReg.dPhi = hJ0.DeltaPhi(hJ1)
405     HNoReg.dEta = abs(hJ0.Eta()-hJ1.Eta())
406 nmohr 1.21 hJet_MtArray[0][0] = hJ0.Mt()
407     hJet_MtArray[1][0] = hJ1.Mt()
408 nmohr 1.24 hJet_EtArray[0][0] = hJ0.Et()
409     hJet_EtArray[1][0] = hJ1.Et()
410 peller 1.15 rPt0 = max(0.0001,readerJet0.EvaluateRegression( "jet0Regression" )[0])
411     rPt1 = max(0.0001,readerJet1.EvaluateRegression( "jet1Regression" )[0])
412 nmohr 1.1 hJet_regWeight[0] = rPt0/hJet_pt0
413     hJet_regWeight[1] = rPt1/hJet_pt1
414     rE0 = hJet_e0*hJet_regWeight[0]
415     rE1 = hJet_e1*hJet_regWeight[1]
416     hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
417     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
418     tree.hJet_pt[0] = rPt0
419     tree.hJet_pt[1] = rPt1
420     tree.hJet_e[0] = rE0
421     tree.hJet_e[1] = rE1
422     H.HiggsFlag = 1
423     H.mass = (hJ0+hJ1).M()
424     H.pt = (hJ0+hJ1).Pt()
425     H.eta = (hJ0+hJ1).Eta()
426     H.phi = (hJ0+hJ1).Phi()
427     H.dR = hJ0.DeltaR(hJ1)
428     H.dPhi = hJ0.DeltaPhi(hJ1)
429     H.dEta = abs(hJ0.Eta()-hJ1.Eta())
430     if hJet_regWeight[0] > 5. or hJet_regWeight[1] > 5.:
431 nmohr 1.24 print 'Event %.0f' %(Event[0])
432 nmohr 1.1 print 'MET %.2f' %(METet[0])
433     print 'rho25 %.2f' %(rho25[0])
434     for key, value in regDict.items():
435 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'):
436 nmohr 1.24 print '%s 0: %.2f'%(key, theVars0[key][0])
437 nmohr 1.26 print '%s 1: %.2f'%(key, theVars1[key][0])
438 nmohr 1.1 for i in range(2):
439     print 'dPhi %.0f %.2f' %(i,hJet_MET_dPhiArray[i][0])
440 nmohr 1.24 for i in range(2):
441     print 'ptRaw %.0f %.2f' %(i,hJet_ptRawArray[i][0])
442     for i in range(2):
443     print 'Mt %.0f %.2f' %(i,hJet_MtArray[i][0])
444     for i in range(2):
445     print 'Et %.0f %.2f' %(i,hJet_EtArray[i][0])
446 nmohr 1.1 print 'corr 0 %.2f' %(hJet_regWeight[0])
447     print 'corr 1 %.2f' %(hJet_regWeight[1])
448     print 'rPt0 %.2f' %(rPt0)
449     print 'rPt1 %.2f' %(rPt1)
450     print 'rE0 %.2f' %(rE0)
451     print 'rE1 %.2f' %(rE1)
452     print 'Mass %.2f' %(H.mass)
453 nmohr 1.34 if fatHiggsFlag:
454     hFJ0.SetPtEtaPhiE(fathFilterJets_pt0,fathFilterJets_eta0,fathFilterJets_phi0,fathFilterJets_e0)
455     hFJ1.SetPtEtaPhiE(fathFilterJets_pt1,fathFilterJets_eta1,fathFilterJets_phi1,fathFilterJets_e1)
456     rFJPt0 = max(0.0001,readerFJ0.EvaluateRegression( "jet0RegressionFJ" )[0])
457     rFJPt1 = max(0.0001,readerFJ1.EvaluateRegression( "jet1RegressionFJ" )[0])
458     fathFilterJets_regWeight[0] = rPt0/fathFilterJets_pt0
459     fathFilterJets_regWeight[1] = rPt1/fathFilterJets_pt1
460     rFJE0 = fathFilterJets_e0*fathFilterJets_regWeight[0]
461     rFJE1 = fathFilterJets_e1*fathFilterJets_regWeight[1]
462     hFJ0.SetPtEtaPhiE(rFJPt0,fathFilterJets_eta0,fathFilterJets_phi0,rFJE0)
463     hFJ1.SetPtEtaPhiE(rFJPt1,fathFilterJets_eta1,fathFilterJets_phi1,rFJE1)
464     FatHReg[0] = (hFJ0+hFJ1).M()
465     FatHReg[1] = (hFJ0+hFJ1).Pt()
466     else:
467     FatHReg[0] = 0.
468     FatHReg[1] = 0.
469    
470     #print rFJPt0
471     #print rFJPt1
472 nmohr 1.30
473     angleHB[0]=fAngleHB.EvalInstance()
474     angleLZ[0]=fAngleLZ.EvalInstance()
475     angleZZS[0]=fAngleZZS.EvalInstance()
476    
477     for i, angLikeBkg in enumerate(AngLikeBkgs):
478     likeSBH[i] = math.fabs(SigBH[i].Eval(angleHB[0]))
479     likeBBH[i] = math.fabs(BkgBH[i].Eval(angleHB[0]))
480    
481     likeSZZS[i] = math.fabs(SigZZS[i].Eval(angleZZS[0]))
482     likeBZZS[i] = math.fabs(BkgZZS[i].Eval(angleZZS[0]))
483    
484     likeSLZ[i] = math.fabs(SigLZ[i].Eval(angleLZ[0]))
485     likeBLZ[i] = math.fabs(BkgLZ[i].Eval(angleLZ[0]))
486    
487     likeSMassZS[i] = math.fabs(SigMassZS[i].Eval(fHVMass.EvalInstance()))
488     likeBMassZS[i] = math.fabs(BkgMassZS[i].Eval(fHVMass.EvalInstance()))
489    
490     scaleSig = float( ang_yield['Sig'] / (ang_yield['Sig'] + ang_yield[angLikeBkg]))
491     scaleBkg = float( ang_yield[angLikeBkg] / (ang_yield['Sig'] + ang_yield[angLikeBkg]) )
492    
493     numerator = (likeSBH[i]*likeSZZS[i]*likeSLZ[i]*likeSMassZS[i]);
494     denominator = ((scaleBkg*likeBLZ[i]*likeBZZS[i]*likeBBH[i]*likeBMassZS[i])+(scaleSig*likeSBH[i]*likeSZZS[i]*likeSLZ[i]*likeSMassZS[i]))
495    
496     if denominator > 0:
497     kinLikeRatio[i] = numerator/denominator;
498     else:
499     kinLikeRatio[i] = 0;
500    
501 nmohr 1.1 if job.type == 'DATA':
502     newtree.Fill()
503     continue
504    
505     for i in range(2):
506 peller 1.15 flavour = int(tree.hJet_flavour[i])
507     pt = float(tree.hJet_pt[i])
508     eta = float(tree.hJet_eta[i])
509     csv = float(tree.hJet_csv[i])
510 nmohr 1.26 hJet_csvOld[i] = csv
511     if anaTag == '7TeV':
512     tree.hJet_csv[i] = corrCSV(btagNom,csv,flavour)
513     hJet_csvDown[i] = corrCSV(btagDown,csv,flavour)
514     hJet_csvUp[i] = corrCSV(btagUp,csv,flavour)
515     hJet_csvFDown[i] = corrCSV(btagFDown,csv,flavour)
516     hJet_csvFUp[i] = corrCSV(btagFUp,csv,flavour)
517     elif anaTag == '8TeV':
518     tree.hJet_csv[i] = btagNom.reshape(eta,pt,csv,flavour)
519     hJet_csvDown[i] = btagDown.reshape(eta,pt,csv,flavour)
520     hJet_csvUp[i] = btagUp.reshape(eta,pt,csv,flavour)
521     hJet_csvFDown[i] = btagFDown.reshape(eta,pt,csv,flavour)
522     hJet_csvFUp[i] = btagFUp.reshape(eta,pt,csv,flavour)
523 nmohr 1.1
524     for updown in ['up','down']:
525     #JER
526     if updown == 'up':
527     inner = 0.06
528     outer = 0.1
529     if updown == 'down':
530     inner = -0.06
531     outer = -0.1
532     #Calculate
533     if abs(hJet_eta0)<1.1: res0 = inner
534     else: res0 = outer
535     if abs(hJet_eta1)<1.1: res1 = inner
536     else: res1 = outer
537     rPt0 = hJet_pt0 + (hJet_pt0-hJet_genPt0)*res0
538     rPt1 = hJet_pt1 + (hJet_pt1-hJet_genPt1)*res1
539     rE0 = hJet_e0*rPt0/hJet_pt0
540     rE1 = hJet_e1*rPt1/hJet_pt1
541 nmohr 1.27 if applyRegression:
542     for key in regVars:
543     var = regDict[key]
544     if var == 'Jet_pt' or var == 'Jet_e' or var == 'hJet_pt' or var == 'hJet_e' or var == 'Jet_ptRaw':
545     if var == 'Jet_ptRaw':
546     hJet_ptRawArray[0][0] = hJet_ptRaw0*corrRes0*rPt0/hJet_pt0
547     hJet_ptRawArray[1][0] = hJet_ptRaw1*rPt1/hJet_pt1
548    
549     elif var == 'Jet_pt' or var == 'hJet_pt':
550     theVars0[key] = rPt0
551     theVars1[key] = rPt1
552     elif var == 'Jet_e' or var == 'hJet_e':
553     theVars0[key] = rE0
554     theVars1[key] = rE1
555     hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
556     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
557     hJet_MtArray[0][0] = hJ0.Mt()
558     hJet_MtArray[1][0] = hJ1.Mt()
559     hJet_EtArray[0][0] = hJ0.Et()
560     hJet_EtArray[1][0] = hJ1.Et()
561     rPt0 = max(0.0001,readerJet0.EvaluateRegression( "jet0Regression" )[0])
562     rPt1 = max(0.0001,readerJet1.EvaluateRegression( "jet1Regression" )[0])
563     rE0 = hJet_e0*rPt0/hJet_pt0
564     rE1 = hJet_e1*rPt1/hJet_pt1
565 nmohr 1.1 hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
566     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
567     #Set
568     if updown == 'up':
569     hJet_pt_JER_up[0]=rPt0
570     hJet_pt_JER_up[1]=rPt1
571     hJet_e_JER_up[0]=rE0
572     hJet_e_JER_up[1]=rE1
573     H_JER[0]=(hJ0+hJ1).M()
574     H_JER[2]=(hJ0+hJ1).Pt()
575     if updown == 'down':
576     hJet_pt_JER_down[0]=rPt0
577     hJet_pt_JER_down[1]=rPt1
578     hJet_e_JER_down[0]=rE0
579     hJet_e_JER_down[1]=rE1
580     H_JER[1]=(hJ0+hJ1).M()
581     H_JER[3]=(hJ0+hJ1).Pt()
582    
583     #JES
584     if updown == 'up':
585     variation=1
586     if updown == 'down':
587     variation=-1
588     #calculate
589     rPt0 = hJet_pt0*(1+variation*hJet_JECUnc0)
590     rPt1 = hJet_pt1*(1+variation*hJet_JECUnc1)
591     rE0 = hJet_e0*(1+variation*hJet_JECUnc0)
592     rE1 = hJet_e1*(1+variation*hJet_JECUnc1)
593 nmohr 1.27 if applyRegression:
594     for key in regVars:
595     var = regDict[key]
596     if var == 'Jet_pt' or var == 'Jet_e' or var == 'hJet_pt' or var == 'hJet_e' or var == 'Jet_ptRaw':
597     if var == 'Jet_ptRaw':
598     hJet_ptRawArray[0][0] = hJet_ptRaw0*(1+variation*hJet_JECUnc0)
599     hJet_ptRawArray[1][0] = hJet_ptRaw1*(1+variation*hJet_JECUnc1)
600    
601     elif var == 'Jet_pt' or var == 'hJet_pt':
602     theVars0[key] = rPt0
603     theVars1[key] = rPt1
604     elif var == 'Jet_e' or var == 'hJet_e':
605     theVars0[key] = rE0
606     theVars1[key] = rE1
607     hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
608     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
609     hJet_MtArray[0][0] = hJ0.Mt()
610     hJet_MtArray[1][0] = hJ1.Mt()
611     hJet_EtArray[0][0] = hJ0.Et()
612     hJet_EtArray[1][0] = hJ1.Et()
613     rPt0 = max(0.0001,readerJet0.EvaluateRegression( "jet0Regression" )[0])
614     rPt1 = max(0.0001,readerJet1.EvaluateRegression( "jet1Regression" )[0])
615     rE0 = hJet_e0*rPt0/hJet_pt0
616     rE1 = hJet_e1*rPt1/hJet_pt1
617 nmohr 1.1 hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
618     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
619     #Fill
620     if updown == 'up':
621     hJet_pt_JES_up[0]=rPt0
622     hJet_pt_JES_up[1]=rPt1
623     hJet_e_JES_up[0]=rE0
624     hJet_e_JES_up[1]=rE1
625     H_JES[0]=(hJ0+hJ1).M()
626     H_JES[2]=(hJ0+hJ1).Pt()
627     if updown == 'down':
628     hJet_pt_JES_down[0]=rPt0
629     hJet_pt_JES_down[1]=rPt1
630     hJet_e_JES_down[0]=rE0
631     hJet_e_JES_down[1]=rE1
632     H_JES[1]=(hJ0+hJ1).M()
633     H_JES[3]=(hJ0+hJ1).Pt()
634    
635     newtree.Fill()
636    
637     newtree.AutoSave()
638     output.Close()