ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/write_regression_systematics.py
Revision: 1.38
Committed: Wed Mar 20 11:29:47 2013 UTC (12 years, 1 month ago) by nmohr
Content type: text/x-python
Branch: MAIN
CVS Tags: LHCP_PreAppFreeze
Changes since 1.37: +1 -1 lines
Log Message:
Fix

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 nmohr 1.37 if(eta< 0.5): return 0.052
67     if(eta< 1.1): return 0.057
68     if(eta< 1.7): return 0.096
69     if(eta< 2.3): return 0.134
70     if(eta< 5): return 0.28
71 peller 1.15 return 0
72    
73     def corrPt(pt,eta,mcPt):
74     return (pt+resolutionBias(math.fabs(eta))*(pt-mcPt))/pt
75    
76 nmohr 1.1 def corrCSV(btag, csv, flav):
77     if(csv < 0.): return csv
78     if(csv > 1.): return csv;
79     if(flav == 0): return csv;
80     if(math.fabs(flav) == 5): return btag.ib.Eval(csv)
81     if(math.fabs(flav) == 4): return btag.ic.Eval(csv)
82     if(math.fabs(flav) != 4 and math.fabs(flav) != 5): return btag.il.Eval(csv)
83     return -10000
84    
85    
86 peller 1.15 def csvReshape(sh, pt, eta, csv, flav):
87     return sh.reshape(float(eta), float(pt), float(csv), int(flav))
88    
89    
90 nmohr 1.1 for job in info:
91 peller 1.15 if not job.name in namelist: 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 nmohr 1.36 else:
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 nmohr 1.36 lhe_weight_map = False if not config.has_option('LHEWeights', 'weights_per_bin') else eval(config.get('LHEWeights', 'weights_per_bin'))
127    
128    
129 nmohr 1.1 print '\t - %s' %(job.name)
130 nmohr 1.34 input = ROOT.TFile.Open(pathIN+job.get_path,'read')
131 nmohr 1.35 output = ROOT.TFile.Open(pathOUT+job.get_path,'recreate')
132 nmohr 1.1
133     input.cd()
134 nmohr 1.36 if lhe_weight_map and 'DY' in job.name:
135     inclusiveJob = info.get_sample('DY')
136     print inclusiveJob.name
137     inclusive = ROOT.TFile.Open(pathIN+inclusiveJob.get_path,'read')
138     inclusive.cd()
139     obj = ROOT.TObject
140     for key in ROOT.gDirectory.GetListOfKeys():
141     input.cd()
142     obj = key.ReadObj()
143     if obj.GetName() == job.tree:
144     continue
145     output.cd()
146     obj.Write(key.GetName())
147     inclusive.Close()
148     else:
149     obj = ROOT.TObject
150     for key in ROOT.gDirectory.GetListOfKeys():
151     input.cd()
152     obj = key.ReadObj()
153     if obj.GetName() == job.tree:
154     continue
155     output.cd()
156     obj.Write(key.GetName())
157 nmohr 1.1
158 nmohr 1.3 input.cd()
159 nmohr 1.1 tree = input.Get(job.tree)
160     nEntries = tree.GetEntries()
161    
162     H = ROOT.H()
163     HNoReg = ROOT.H()
164     tree.SetBranchStatus('H',0)
165     output.cd()
166     newtree = tree.CloneTree(0)
167    
168     hJ0 = ROOT.TLorentzVector()
169     hJ1 = ROOT.TLorentzVector()
170 nmohr 1.37 vect = ROOT.TLorentzVector()
171 nmohr 1.36 #hFJ0 = ROOT.TLorentzVector()
172     #hFJ1 = ROOT.TLorentzVector()
173 nmohr 1.1
174 peller 1.15 regWeight = config.get("Regression","regWeight")
175     regDict = eval(config.get("Regression","regDict"))
176     regVars = eval(config.get("Regression","regVars"))
177 nmohr 1.36 #regWeightFilterJets = config.get("Regression","regWeightFilterJets")
178     #regDictFilterJets = eval(config.get("Regression","regDictFilterJets"))
179     #regVarsFilterJets = eval(config.get("Regression","regVarsFilterJets"))
180 nmohr 1.1
181     #Regression branches
182     applyRegression = True
183     hJet_pt = array('f',[0]*2)
184     hJet_e = array('f',[0]*2)
185 nmohr 1.34 newtree.Branch( 'H', H , 'HiggsFlag/I:mass/F:pt/F:eta/F:phi/F:dR/F:dPhi/F:dEta/F' )
186     newtree.Branch( 'HNoReg', HNoReg , 'HiggsFlag/I:mass/F:pt/F:eta/F:phi/F:dR/F:dPhi/F:dEta/F' )
187 nmohr 1.36 #FatHReg = array('f',[0]*2)
188     #newtree.Branch('FatHReg',FatHReg,'filteredmass:filteredpt/F')
189 nmohr 1.1 Event = array('f',[0])
190     METet = array('f',[0])
191     rho25 = array('f',[0])
192     METphi = array('f',[0])
193     fRho25 = ROOT.TTreeFormula("rho25",'rho25',tree)
194     fEvent = ROOT.TTreeFormula("Event",'EVENT.event',tree)
195 nmohr 1.34 fFatHFlag = ROOT.TTreeFormula("FatHFlag",'FatH.FatHiggsFlag',tree)
196     fFatHnFilterJets = ROOT.TTreeFormula("FatHnFilterJets",'nfathFilterJets',tree)
197 nmohr 1.1 fMETet = ROOT.TTreeFormula("METet",'METnoPU.et',tree)
198     fMETphi = ROOT.TTreeFormula("METphi",'METnoPU.phi',tree)
199 nmohr 1.30 fHVMass = ROOT.TTreeFormula("HVMass",'HVMass',tree)
200 nmohr 1.21 hJet_MtArray = [array('f',[0]),array('f',[0])]
201 nmohr 1.24 hJet_EtArray = [array('f',[0]),array('f',[0])]
202 nmohr 1.1 hJet_MET_dPhi = array('f',[0]*2)
203     hJet_regWeight = array('f',[0]*2)
204 nmohr 1.34 fathFilterJets_regWeight = array('f',[0]*2)
205 nmohr 1.1 hJet_MET_dPhiArray = [array('f',[0]),array('f',[0])]
206 nmohr 1.27 hJet_ptRawArray = [array('f',[0]),array('f',[0])]
207 nmohr 1.1 newtree.Branch('hJet_MET_dPhi',hJet_MET_dPhi,'hJet_MET_dPhi[2]/F')
208     newtree.Branch('hJet_regWeight',hJet_regWeight,'hJet_regWeight[2]/F')
209     readerJet0 = ROOT.TMVA.Reader("!Color:!Silent" )
210     readerJet1 = ROOT.TMVA.Reader("!Color:!Silent" )
211    
212     theForms = {}
213     theVars0 = {}
214     theVars1 = {}
215 nmohr 1.34 def addVarsToReader(reader,regDict,regVars,theVars,theForms,i,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray):
216 nmohr 1.24 for key in regVars:
217     var = regDict[key]
218     theVars[key] = array( 'f', [ 0 ] )
219     if var == 'Jet_MET_dPhi':
220     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
221     reader.AddVariable( key, hJet_MET_dPhiArray[i] )
222     elif var == 'METet':
223     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
224 nmohr 1.26 reader.AddVariable( key, METet )
225 nmohr 1.24 elif var == 'rho25':
226     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
227 nmohr 1.26 reader.AddVariable( key, rho25 )
228 nmohr 1.24 elif var == 'Jet_mt':
229     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
230 nmohr 1.26 reader.AddVariable( key, hJet_MtArray[i] )
231 nmohr 1.24 elif var == 'Jet_et':
232     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
233 nmohr 1.26 reader.AddVariable( key, hJet_EtArray[i] )
234 nmohr 1.24 elif var == 'Jet_ptRaw':
235     print 'Adding var: %s with %s to readerJet%.0f' %(key,var,i)
236 nmohr 1.26 reader.AddVariable( key, hJet_ptRawArray[i] )
237 nmohr 1.24 else:
238     reader.AddVariable(key,theVars[key])
239     formula = regDict[key].replace("[0]","[%.0f]" %i)
240     print 'Adding var: %s with %s to readerJet%.0f' %(key,formula,i)
241     theForms['form_reg_%s_%.0f'%(key,i)] = ROOT.TTreeFormula("form_reg_%s_%.0f"%(key,i),'%s' %(formula),tree)
242     return
243 nmohr 1.34 addVarsToReader(readerJet0,regDict,regVars,theVars0,theForms,0,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray)
244     addVarsToReader(readerJet1,regDict,regVars,theVars1,theForms,1,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray)
245 nmohr 1.24 readerJet0.BookMVA( "jet0Regression", regWeight )
246     readerJet1.BookMVA( "jet1Regression", regWeight )
247    
248 nmohr 1.36 #readerFJ0 = ROOT.TMVA.Reader("!Color:!Silent" )
249     #readerFJ1 = ROOT.TMVA.Reader("!Color:!Silent" )
250     #theFormsFJ = {}
251     #theVars0FJ = {}
252     #theVars1FJ = {}
253     #addVarsToReader(readerFJ0,regDictFilterJets,regVarsFilterJets,theVars0FJ,theFormsFJ,0,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray)
254     #addVarsToReader(readerFJ1,regDictFilterJets,regVarsFilterJets,theVars1FJ,theFormsFJ,1,hJet_MET_dPhiArray,METet,rho25,hJet_MtArray,hJet_EtArray,hJet_ptRawArray)
255     #readerFJ0.BookMVA( "jet0RegressionFJ", regWeightFilterJets )
256     #readerFJ1.BookMVA( "jet1RegressionFJ", regWeightFilterJets )
257 nmohr 1.24
258 nmohr 1.1
259 bortigno 1.4 #Add training Flag
260 nmohr 1.18 EventForTraining = array('i',[0])
261     newtree.Branch('EventForTraining',EventForTraining,'EventForTraining/I')
262 bortigno 1.4 EventForTraining[0]=0
263 bortigno 1.6
264 bortigno 1.5 TFlag=ROOT.TTreeFormula("EventForTraining","EVENT.event%2",tree)
265 nmohr 1.30
266     #Angular Likelihood
267     angleHB = array('f',[0])
268     newtree.Branch('angleHB',angleHB,'angleHB/F')
269     angleLZ = array('f',[0])
270     newtree.Branch('angleLZ',angleLZ,'angleLZ/F')
271     angleZZS = array('f',[0])
272     newtree.Branch('angleZZS',angleZZS,'angleZZS/F')
273     kinLikeRatio = array('f',[0]*len(AngLikeBkgs))
274     newtree.Branch('kinLikeRatio',kinLikeRatio,'%s/F' %(':'.join(AngLikeBkgs)))
275 nmohr 1.37 fAngleHB = ROOT.TTreeFormula("fAngleHB",'abs(VHbb::ANGLEHB(hJet_pt[0],hJet_eta[0],hJet_phi[0],hJet_e[0],hJet_pt[1],hJet_eta[1],hJet_phi[1],hJet_e[1]))',newtree)
276 nmohr 1.30 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)
277     fAngleZZS = ROOT.TTreeFormula("fAngleZZS",'abs(VHbb::ANGLELZ(H.pt,H.eta,H.phi,H.pt,V.pt,V.eta,V.phi,V.mass))',newtree)
278 nmohr 1.37 fVpt = ROOT.TTreeFormula("fVpt",'V.pt',tree)
279     fVeta = ROOT.TTreeFormula("fVeta",'V.eta',tree)
280     fVphi = ROOT.TTreeFormula("fVphi",'V.phi',tree)
281     fVmass = ROOT.TTreeFormula("fVmass",'V.mass',tree)
282 nmohr 1.30 likeSBH = array('f',[0]*len(AngLikeBkgs))
283     likeBBH = array('f',[0]*len(AngLikeBkgs))
284     likeSLZ = array('f',[0]*len(AngLikeBkgs))
285     likeBLZ = array('f',[0]*len(AngLikeBkgs))
286     likeSZZS = array('f',[0]*len(AngLikeBkgs))
287     likeBZZS = array('f',[0]*len(AngLikeBkgs))
288     likeSMassZS = array('f',[0]*len(AngLikeBkgs))
289     likeBMassZS = array('f',[0]*len(AngLikeBkgs))
290 nmohr 1.37 HVMass_Reg = array('f',[0])
291     newtree.Branch('HVMass_Reg',HVMass_Reg,'HVMass_Reg/F')
292 nmohr 1.30
293     SigBH = []; BkgBH = []; SigLZ = []; BkgLZ = []; SigZZS = []; BkgZZS = []; SigMassZS = []; BkgMassZS = [];
294     for angLikeBkg in AngLikeBkgs:
295     f = ROOT.TFile("../data/angleFitFunctions-%s.root"%(angLikeBkg),"READ")
296     SigBH.append(f.Get("sigFuncBH"))
297     BkgBH.append(f.Get("bkgFuncBH"))
298     SigLZ.append(f.Get("sigFuncLZ"))
299     BkgLZ.append(f.Get("bkgFuncLZ"))
300     SigZZS.append(f.Get("sigFuncZZS"))
301     BkgZZS.append(f.Get("bkgFuncZZS"))
302     SigMassZS.append(f.Get("sigFuncMassZS"))
303     BkgMassZS.append(f.Get("bkgFuncMassZS"))
304     f.Close()
305    
306 nmohr 1.1
307     if job.type != 'DATA':
308     #CSV branches
309     hJet_flavour = array('f',[0]*2)
310     hJet_csv = array('f',[0]*2)
311     hJet_csvOld = array('f',[0]*2)
312     hJet_csvUp = array('f',[0]*2)
313     hJet_csvDown = array('f',[0]*2)
314     hJet_csvFUp = array('f',[0]*2)
315     hJet_csvFDown = array('f',[0]*2)
316     newtree.Branch('hJet_csvOld',hJet_csvOld,'hJet_csvOld[2]/F')
317     newtree.Branch('hJet_csvUp',hJet_csvUp,'hJet_csvUp[2]/F')
318     newtree.Branch('hJet_csvDown',hJet_csvDown,'hJet_csvDown[2]/F')
319     newtree.Branch('hJet_csvFUp',hJet_csvFUp,'hJet_csvFUp[2]/F')
320     newtree.Branch('hJet_csvFDown',hJet_csvFDown,'hJet_csvFDown[2]/F')
321    
322     #JER branches
323     hJet_pt_JER_up = array('f',[0]*2)
324     newtree.Branch('hJet_pt_JER_up',hJet_pt_JER_up,'hJet_pt_JER_up[2]/F')
325     hJet_pt_JER_down = array('f',[0]*2)
326     newtree.Branch('hJet_pt_JER_down',hJet_pt_JER_down,'hJet_pt_JER_down[2]/F')
327     hJet_e_JER_up = array('f',[0]*2)
328     newtree.Branch('hJet_e_JER_up',hJet_e_JER_up,'hJet_e_JER_up[2]/F')
329     hJet_e_JER_down = array('f',[0]*2)
330     newtree.Branch('hJet_e_JER_down',hJet_e_JER_down,'hJet_e_JER_down[2]/F')
331     H_JER = array('f',[0]*4)
332     newtree.Branch('H_JER',H_JER,'mass_up:mass_down:pt_up:pt_down/F')
333 nmohr 1.37 HVMass_JER_up = array('f',[0])
334     HVMass_JER_down = array('f',[0])
335     newtree.Branch('HVMass_JER_up',HVMass_JER_up,'HVMass_JER_up/F')
336     newtree.Branch('HVMass_JER_down',HVMass_JER_down,'HVMass_JER_down/F')
337     angleHB_JER_up = array('f',[0])
338     angleHB_JER_down = array('f',[0])
339     angleZZS_JER_up = array('f',[0])
340     angleZZS_JER_down = array('f',[0])
341     newtree.Branch('angleHB_JER_up',angleHB_JER_up,'angleHB_JER_up/F')
342     newtree.Branch('angleHB_JER_down',angleHB_JER_down,'angleHB_JER_down/F')
343     newtree.Branch('angleZZS_JER_up',angleZZS_JER_up,'angleZZS_JER_up/F')
344     newtree.Branch('angleZZS_JER_down',angleZZS_JER_down,'angleZZS_JER_down/F')
345 nmohr 1.1
346     #JES branches
347     hJet_pt_JES_up = array('f',[0]*2)
348     newtree.Branch('hJet_pt_JES_up',hJet_pt_JES_up,'hJet_pt_JES_up[2]/F')
349     hJet_pt_JES_down = array('f',[0]*2)
350     newtree.Branch('hJet_pt_JES_down',hJet_pt_JES_down,'hJet_pt_JES_down[2]/F')
351     hJet_e_JES_up = array('f',[0]*2)
352     newtree.Branch('hJet_e_JES_up',hJet_e_JES_up,'hJet_e_JES_up[2]/F')
353     hJet_e_JES_down = array('f',[0]*2)
354     newtree.Branch('hJet_e_JES_down',hJet_e_JES_down,'hJet_e_JES_down[2]/F')
355     H_JES = array('f',[0]*4)
356     newtree.Branch('H_JES',H_JES,'mass_up:mass_down:pt_up:pt_down/F')
357 nmohr 1.37 HVMass_JES_up = array('f',[0])
358     HVMass_JES_down = array('f',[0])
359     newtree.Branch('HVMass_JES_up',HVMass_JES_up,'HVMass_JES_up/F')
360     newtree.Branch('HVMass_JES_down',HVMass_JES_down,'HVMass_JES_down/F')
361     angleHB_JES_up = array('f',[0])
362     angleHB_JES_down = array('f',[0])
363     angleZZS_JES_up = array('f',[0])
364     angleZZS_JES_down = array('f',[0])
365     newtree.Branch('angleHB_JES_up',angleHB_JES_up,'angleHB_JES_up/F')
366     newtree.Branch('angleHB_JES_down',angleHB_JES_down,'angleHB_JES_down/F')
367     newtree.Branch('angleZZS_JES_up',angleZZS_JES_up,'angleZZS_JES_up/F')
368     newtree.Branch('angleZZS_JES_down',angleZZS_JES_down,'angleZZS_JES_down/F')
369    
370     #Formulas for syst in angular
371     fAngleHB_JER_up = ROOT.TTreeFormula("fAngleHB_JER_up",'abs(VHbb::ANGLEHB(hJet_pt_JER_up[0],hJet_eta[0],hJet_phi[0],hJet_e_JER_up[0],hJet_pt_JER_up[1],hJet_eta[1],hJet_phi[1],hJet_e_JER_up[1]))',newtree)
372     fAngleHB_JER_down = ROOT.TTreeFormula("fAngleHB_JER_down",'abs(VHbb::ANGLEHB(hJet_pt_JER_down[0],hJet_eta[0],hJet_phi[0],hJet_e_JER_down[0],hJet_pt_JER_down[1],hJet_eta[1],hJet_phi[1],hJet_e_JER_down[1]))',newtree)
373     fAngleHB_JES_up = ROOT.TTreeFormula("fAngleHB_JES_up",'abs(VHbb::ANGLEHB(hJet_pt_JES_up[0],hJet_eta[0],hJet_phi[0],hJet_e_JES_up[0],hJet_pt_JES_up[1],hJet_eta[1],hJet_phi[1],hJet_e_JES_up[1]))',newtree)
374     fAngleHB_JES_down = ROOT.TTreeFormula("fAngleHB_JES_down",'abs(VHbb::ANGLEHB(hJet_pt_JES_down[0],hJet_eta[0],hJet_phi[0],hJet_e_JES_down[0],hJet_pt_JES_down[1],hJet_eta[1],hJet_phi[1],hJet_e_JES_down[1]))',newtree)
375     fAngleZZS_JER_up = ROOT.TTreeFormula("fAngleZZS_JER_Up",'abs(VHbb::ANGLELZ(H_JER.pt_up,H.eta,H.phi,H_JER.pt_up,V.pt,V.eta,V.phi,V.mass))',newtree)
376     fAngleZZS_JER_down = ROOT.TTreeFormula("fAngleZZS_JER_Down",'abs(VHbb::ANGLELZ(H_JER.pt_down,H.eta,H.phi,H_JER.pt_down,V.pt,V.eta,V.phi,V.mass))',newtree)
377     fAngleZZS_JES_up = ROOT.TTreeFormula("fAngleZZS_JES_Up",'abs(VHbb::ANGLELZ(H_JER.pt_up,H.eta,H.phi,H_JER.pt_up,V.pt,V.eta,V.phi,V.mass))',newtree)
378     fAngleZZS_JES_down = ROOT.TTreeFormula("fAngleZZS_JES_Down",'abs(VHbb::ANGLELZ(H_JER.pt_down,H.eta,H.phi,H_JER.pt_down,V.pt,V.eta,V.phi,V.mass))',newtree)
379 nmohr 1.27 lheWeight = array('f',[0])
380 nmohr 1.36 newtree.Branch('lheWeight',lheWeight,'lheWeight/F')
381     theBinForms = {}
382     if lhe_weight_map and 'DY' in job.name:
383     for bin in lhe_weight_map:
384     theBinForms[bin] = ROOT.TTreeFormula("Bin_formula_%s"%(bin),bin,tree)
385     else:
386 nmohr 1.26 lheWeight[0] = 1.
387 nmohr 1.1
388     #iter=0
389    
390    
391     for entry in range(0,nEntries):
392     tree.GetEntry(entry)
393    
394 nmohr 1.38 if job.type != 'DATA':
395 nmohr 1.18 EventForTraining[0]=int(not TFlag.EvalInstance())
396 nmohr 1.36 if lhe_weight_map and 'DY' in job.name:
397     match_bin = None
398     for bin in lhe_weight_map:
399     if theBinForms[bin].EvalInstance() > 0.:
400     match_bin = bin
401     if match_bin:
402     lheWeight[0] = lhe_weight_map[match_bin]
403     else:
404     lheWeight[0] = 1.
405 bortigno 1.9
406 nmohr 1.34 #Has fat higgs
407 nmohr 1.36 #fatHiggsFlag=fFatHFlag.EvalInstance()*fFatHnFilterJets.EvalInstance()
408 nmohr 1.1
409     #get
410 nmohr 1.37 vect.SetPtEtaPhiM(fVpt.EvalInstance(),fVeta.EvalInstance(),fVphi.EvalInstance(),fVmass.EvalInstance())
411 nmohr 1.1 hJet_pt = tree.hJet_pt
412     hJet_e = tree.hJet_e
413     hJet_pt0 = tree.hJet_pt[0]
414     hJet_pt1 = tree.hJet_pt[1]
415     hJet_eta0 = tree.hJet_eta[0]
416     hJet_eta1 = tree.hJet_eta[1]
417     hJet_genPt0 = tree.hJet_genPt[0]
418     hJet_genPt1 = tree.hJet_genPt[1]
419 nmohr 1.27 hJet_ptRaw0 = tree.hJet_ptRaw[0]
420     hJet_ptRaw1 = tree.hJet_ptRaw[1]
421 nmohr 1.1 hJet_e0 = tree.hJet_e[0]
422     hJet_e1 = tree.hJet_e[1]
423     hJet_phi0 = tree.hJet_phi[0]
424     hJet_phi1 = tree.hJet_phi[1]
425     hJet_JECUnc0 = tree.hJet_JECUnc[0]
426     hJet_JECUnc1 = tree.hJet_JECUnc[1]
427 nmohr 1.34 #Filterjets
428 nmohr 1.36 #if fatHiggsFlag:
429     # fathFilterJets_pt0 = tree.fathFilterJets_pt[0]
430     # fathFilterJets_pt1 = tree.fathFilterJets_pt[1]
431     # fathFilterJets_eta0 = tree.fathFilterJets_eta[0]
432     # fathFilterJets_eta1 = tree.fathFilterJets_eta[1]
433     # fathFilterJets_phi0 = tree.fathFilterJets_phi[0]
434     # fathFilterJets_phi1 = tree.fathFilterJets_phi[1]
435     # fathFilterJets_e0 = tree.fathFilterJets_e[0]
436     # fathFilterJets_e1 = tree.fathFilterJets_e[1]
437 nmohr 1.1
438     Event[0]=fEvent.EvalInstance()
439     METet[0]=fMETet.EvalInstance()
440     rho25[0]=fRho25.EvalInstance()
441     METphi[0]=fMETphi.EvalInstance()
442     for key, value in regDict.items():
443 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'):
444 nmohr 1.24 theVars0[key][0] = theForms["form_reg_%s_0" %(key)].EvalInstance()
445     theVars1[key][0] = theForms["form_reg_%s_1" %(key)].EvalInstance()
446 nmohr 1.36 #for key, value in regDictFilterJets.items():
447     # if not (value == 'Jet_MET_dPhi' or value == 'METet' or value == "rho25" or value == "Jet_et" or value == 'Jet_mt' or value == 'Jet_ptRaw'):
448     # theVars0FJ[key][0] = theFormsFJ["form_reg_%s_0" %(key)].EvalInstance()
449     # theVars1FJ[key][0] = theFormsFJ["form_reg_%s_1" %(key)].EvalInstance()
450 nmohr 1.27 hJet_MET_dPhi[0] = deltaPhi(METphi[0],hJet_phi0)
451     hJet_MET_dPhi[1] = deltaPhi(METphi[0],hJet_phi1)
452     hJet_MET_dPhiArray[0][0] = deltaPhi(METphi[0],hJet_phi0)
453     hJet_MET_dPhiArray[1][0] = deltaPhi(METphi[0],hJet_phi1)
454     if not job.type == 'DATA':
455     corrRes0 = corrPt(hJet_pt0,hJet_eta0,hJet_genPt0)
456     corrRes1 = corrPt(hJet_pt1,hJet_eta1,hJet_genPt1)
457     hJet_ptRaw0 *= corrRes0
458     hJet_ptRaw1 *= corrRes1
459     hJet_ptRawArray[0][0] = hJet_ptRaw0
460     hJet_ptRawArray[1][0] = hJet_ptRaw1
461 nmohr 1.37 hJ0.SetPtEtaPhiE(hJet_pt0,hJet_eta0,hJet_phi0,hJet_e0)
462     hJ1.SetPtEtaPhiE(hJet_pt1,hJet_eta1,hJet_phi1,hJet_e1)
463     hJet_et0 = hJ0.Et()
464     hJet_et1 = hJ1.Et()
465     hJet_mt0 = hJ0.Mt()
466     hJet_mt1 = hJ1.Mt()
467 nmohr 1.1
468 peller 1.15
469 nmohr 1.1 if applyRegression:
470     HNoReg.HiggsFlag = 1
471     HNoReg.mass = (hJ0+hJ1).M()
472     HNoReg.pt = (hJ0+hJ1).Pt()
473     HNoReg.eta = (hJ0+hJ1).Eta()
474     HNoReg.phi = (hJ0+hJ1).Phi()
475     HNoReg.dR = hJ0.DeltaR(hJ1)
476     HNoReg.dPhi = hJ0.DeltaPhi(hJ1)
477     HNoReg.dEta = abs(hJ0.Eta()-hJ1.Eta())
478 nmohr 1.21 hJet_MtArray[0][0] = hJ0.Mt()
479     hJet_MtArray[1][0] = hJ1.Mt()
480 nmohr 1.24 hJet_EtArray[0][0] = hJ0.Et()
481     hJet_EtArray[1][0] = hJ1.Et()
482 peller 1.15 rPt0 = max(0.0001,readerJet0.EvaluateRegression( "jet0Regression" )[0])
483     rPt1 = max(0.0001,readerJet1.EvaluateRegression( "jet1Regression" )[0])
484 nmohr 1.1 hJet_regWeight[0] = rPt0/hJet_pt0
485     hJet_regWeight[1] = rPt1/hJet_pt1
486     rE0 = hJet_e0*hJet_regWeight[0]
487     rE1 = hJet_e1*hJet_regWeight[1]
488     hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
489     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
490 nmohr 1.37 #print '###new####'
491     #print 'First regression %s' %rPt0
492 nmohr 1.1 tree.hJet_pt[0] = rPt0
493     tree.hJet_pt[1] = rPt1
494     tree.hJet_e[0] = rE0
495     tree.hJet_e[1] = rE1
496     H.HiggsFlag = 1
497     H.mass = (hJ0+hJ1).M()
498     H.pt = (hJ0+hJ1).Pt()
499     H.eta = (hJ0+hJ1).Eta()
500     H.phi = (hJ0+hJ1).Phi()
501     H.dR = hJ0.DeltaR(hJ1)
502     H.dPhi = hJ0.DeltaPhi(hJ1)
503     H.dEta = abs(hJ0.Eta()-hJ1.Eta())
504 nmohr 1.37 HVMass_Reg[0] = (hJ0+hJ1+vect).M()
505 nmohr 1.1 if hJet_regWeight[0] > 5. or hJet_regWeight[1] > 5.:
506 nmohr 1.24 print 'Event %.0f' %(Event[0])
507 nmohr 1.1 print 'MET %.2f' %(METet[0])
508     print 'rho25 %.2f' %(rho25[0])
509     for key, value in regDict.items():
510 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'):
511 nmohr 1.24 print '%s 0: %.2f'%(key, theVars0[key][0])
512 nmohr 1.26 print '%s 1: %.2f'%(key, theVars1[key][0])
513 nmohr 1.1 for i in range(2):
514     print 'dPhi %.0f %.2f' %(i,hJet_MET_dPhiArray[i][0])
515 nmohr 1.24 for i in range(2):
516     print 'ptRaw %.0f %.2f' %(i,hJet_ptRawArray[i][0])
517     for i in range(2):
518     print 'Mt %.0f %.2f' %(i,hJet_MtArray[i][0])
519     for i in range(2):
520     print 'Et %.0f %.2f' %(i,hJet_EtArray[i][0])
521 nmohr 1.1 print 'corr 0 %.2f' %(hJet_regWeight[0])
522     print 'corr 1 %.2f' %(hJet_regWeight[1])
523     print 'rPt0 %.2f' %(rPt0)
524     print 'rPt1 %.2f' %(rPt1)
525     print 'rE0 %.2f' %(rE0)
526     print 'rE1 %.2f' %(rE1)
527     print 'Mass %.2f' %(H.mass)
528 nmohr 1.36 #if fatHiggsFlag:
529     #hFJ0.SetPtEtaPhiE(fathFilterJets_pt0,fathFilterJets_eta0,fathFilterJets_phi0,fathFilterJets_e0)
530     #hFJ1.SetPtEtaPhiE(fathFilterJets_pt1,fathFilterJets_eta1,fathFilterJets_phi1,fathFilterJets_e1)
531     #rFJPt0 = max(0.0001,readerFJ0.EvaluateRegression( "jet0RegressionFJ" )[0])
532     #rFJPt1 = max(0.0001,readerFJ1.EvaluateRegression( "jet1RegressionFJ" )[0])
533     #fathFilterJets_regWeight[0] = rPt0/fathFilterJets_pt0
534     #fathFilterJets_regWeight[1] = rPt1/fathFilterJets_pt1
535     #rFJE0 = fathFilterJets_e0*fathFilterJets_regWeight[0]
536     #rFJE1 = fathFilterJets_e1*fathFilterJets_regWeight[1]
537     #hFJ0.SetPtEtaPhiE(rFJPt0,fathFilterJets_eta0,fathFilterJets_phi0,rFJE0)
538     #hFJ1.SetPtEtaPhiE(rFJPt1,fathFilterJets_eta1,fathFilterJets_phi1,rFJE1)
539     #FatHReg[0] = (hFJ0+hFJ1).M()
540     #FatHReg[1] = (hFJ0+hFJ1).Pt()
541     #else:
542     #FatHReg[0] = 0.
543     #FatHReg[1] = 0.
544 nmohr 1.34
545     #print rFJPt0
546     #print rFJPt1
547 nmohr 1.30
548     angleHB[0]=fAngleHB.EvalInstance()
549     angleLZ[0]=fAngleLZ.EvalInstance()
550     angleZZS[0]=fAngleZZS.EvalInstance()
551    
552     for i, angLikeBkg in enumerate(AngLikeBkgs):
553     likeSBH[i] = math.fabs(SigBH[i].Eval(angleHB[0]))
554     likeBBH[i] = math.fabs(BkgBH[i].Eval(angleHB[0]))
555    
556     likeSZZS[i] = math.fabs(SigZZS[i].Eval(angleZZS[0]))
557     likeBZZS[i] = math.fabs(BkgZZS[i].Eval(angleZZS[0]))
558    
559     likeSLZ[i] = math.fabs(SigLZ[i].Eval(angleLZ[0]))
560     likeBLZ[i] = math.fabs(BkgLZ[i].Eval(angleLZ[0]))
561    
562     likeSMassZS[i] = math.fabs(SigMassZS[i].Eval(fHVMass.EvalInstance()))
563     likeBMassZS[i] = math.fabs(BkgMassZS[i].Eval(fHVMass.EvalInstance()))
564    
565     scaleSig = float( ang_yield['Sig'] / (ang_yield['Sig'] + ang_yield[angLikeBkg]))
566     scaleBkg = float( ang_yield[angLikeBkg] / (ang_yield['Sig'] + ang_yield[angLikeBkg]) )
567    
568     numerator = (likeSBH[i]*likeSZZS[i]*likeSLZ[i]*likeSMassZS[i]);
569     denominator = ((scaleBkg*likeBLZ[i]*likeBZZS[i]*likeBBH[i]*likeBMassZS[i])+(scaleSig*likeSBH[i]*likeSZZS[i]*likeSLZ[i]*likeSMassZS[i]))
570    
571     if denominator > 0:
572     kinLikeRatio[i] = numerator/denominator;
573     else:
574     kinLikeRatio[i] = 0;
575    
576 nmohr 1.1 if job.type == 'DATA':
577     newtree.Fill()
578     continue
579    
580     for i in range(2):
581 peller 1.15 flavour = int(tree.hJet_flavour[i])
582     pt = float(tree.hJet_pt[i])
583     eta = float(tree.hJet_eta[i])
584     csv = float(tree.hJet_csv[i])
585 nmohr 1.26 hJet_csvOld[i] = csv
586     if anaTag == '7TeV':
587     tree.hJet_csv[i] = corrCSV(btagNom,csv,flavour)
588     hJet_csvDown[i] = corrCSV(btagDown,csv,flavour)
589     hJet_csvUp[i] = corrCSV(btagUp,csv,flavour)
590     hJet_csvFDown[i] = corrCSV(btagFDown,csv,flavour)
591     hJet_csvFUp[i] = corrCSV(btagFUp,csv,flavour)
592 nmohr 1.37 else:
593     #tree.hJet_csv[i] = btagNom.reshape(eta,pt,csv,flavour)
594     #hJet_csvDown[i] = btagDown.reshape(eta,pt,csv,flavour)
595     #hJet_csvUp[i] = btagUp.reshape(eta,pt,csv,flavour)
596     #hJet_csvFDown[i] = btagFDown.reshape(eta,pt,csv,flavour)
597     #hJet_csvFUp[i] = btagFUp.reshape(eta,pt,csv,flavour)
598     tree.hJet_csv[i] = tree.hJet_csv_nominal[i]
599     hJet_csvDown[i] = tree.hJet_csv_downBC[i]
600     hJet_csvUp[i] = tree.hJet_csv_upBC[i]
601     hJet_csvFDown[i] = tree.hJet_csv_downL[i]
602     hJet_csvFUp[i] = tree.hJet_csv_upL[i]
603 nmohr 1.1
604     for updown in ['up','down']:
605     #JER
606     if updown == 'up':
607     inner = 0.06
608     outer = 0.1
609     if updown == 'down':
610     inner = -0.06
611     outer = -0.1
612     #Calculate
613     if abs(hJet_eta0)<1.1: res0 = inner
614     else: res0 = outer
615     if abs(hJet_eta1)<1.1: res1 = inner
616     else: res1 = outer
617     rPt0 = hJet_pt0 + (hJet_pt0-hJet_genPt0)*res0
618     rPt1 = hJet_pt1 + (hJet_pt1-hJet_genPt1)*res1
619     rE0 = hJet_e0*rPt0/hJet_pt0
620     rE1 = hJet_e1*rPt1/hJet_pt1
621 nmohr 1.27 if applyRegression:
622     hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
623     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
624     hJet_MtArray[0][0] = hJ0.Mt()
625     hJet_MtArray[1][0] = hJ1.Mt()
626     hJet_EtArray[0][0] = hJ0.Et()
627     hJet_EtArray[1][0] = hJ1.Et()
628 nmohr 1.37 for key in regVars:
629     var = regDict[key]
630     if key == 'Jet_pt' or key == 'Jet_e' or key == 'hJet_pt' or key == 'hJet_e' or key == 'Jet_ptRaw' or key =='VHbb::evalEt(hJet_pt,hJet_eta,hJet_phi,hJet_e)' or key =='VHbb::evalMt(hJet_pt,hJet_eta,hJet_phi,hJet_e)' or key == 'VHbb::evalJERBias(hJet_ptRaw,hJet_genPt,hJet_eta)':
631     if key == 'Jet_ptRaw':
632     hJet_ptRawArray[0][0] = hJet_ptRaw0*corrRes0*rPt0/hJet_pt0
633     hJet_ptRawArray[1][0] = hJet_ptRaw1*rPt1/hJet_pt1
634     elif key == 'VHbb::evalJERBias(hJet_ptRaw,hJet_genPt,hJet_eta)':
635     theVars0[key][0] = hJet_ptRaw0*corrRes0*rPt0/hJet_pt0
636     theVars1[key][0] = hJet_ptRaw1*rPt1/hJet_pt1
637     elif key == 'Jet_pt' or key == 'hJet_pt':
638     theVars0[key][0] = rPt0
639     theVars1[key][0] = rPt1
640     elif key == 'Jet_e' or key == 'hJet_e':
641     theVars0[key][0] = rE0
642     theVars1[key][0] = rE1
643     elif key == 'VHbb::evalEt(hJet_pt,hJet_eta,hJet_phi,hJet_e)':
644     theVars0[key][0] = hJ0.Et()
645     theVars1[key][0] = hJ1.Et()
646     elif key == 'VHbb::evalMt(hJet_pt,hJet_eta,hJet_phi,hJet_e)':
647     theVars0[key][0] = hJ0.Mt()
648     theVars1[key][0] = hJ1.Mt()
649 nmohr 1.27 rPt0 = max(0.0001,readerJet0.EvaluateRegression( "jet0Regression" )[0])
650     rPt1 = max(0.0001,readerJet1.EvaluateRegression( "jet1Regression" )[0])
651     rE0 = hJet_e0*rPt0/hJet_pt0
652     rE1 = hJet_e1*rPt1/hJet_pt1
653 nmohr 1.1 hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
654     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
655     #Set
656     if updown == 'up':
657     hJet_pt_JER_up[0]=rPt0
658     hJet_pt_JER_up[1]=rPt1
659     hJet_e_JER_up[0]=rE0
660     hJet_e_JER_up[1]=rE1
661     H_JER[0]=(hJ0+hJ1).M()
662     H_JER[2]=(hJ0+hJ1).Pt()
663 nmohr 1.37 HVMass_JER_up[0] = (hJ0+hJ1+vect).M()
664 nmohr 1.1 if updown == 'down':
665     hJet_pt_JER_down[0]=rPt0
666     hJet_pt_JER_down[1]=rPt1
667     hJet_e_JER_down[0]=rE0
668     hJet_e_JER_down[1]=rE1
669     H_JER[1]=(hJ0+hJ1).M()
670     H_JER[3]=(hJ0+hJ1).Pt()
671 nmohr 1.37 HVMass_JER_down[0] = (hJ0+hJ1+vect).M()
672 nmohr 1.1
673     #JES
674     if updown == 'up':
675     variation=1
676     if updown == 'down':
677     variation=-1
678     #calculate
679     rPt0 = hJet_pt0*(1+variation*hJet_JECUnc0)
680     rPt1 = hJet_pt1*(1+variation*hJet_JECUnc1)
681     rE0 = hJet_e0*(1+variation*hJet_JECUnc0)
682     rE1 = hJet_e1*(1+variation*hJet_JECUnc1)
683 nmohr 1.37 #print 'res %s: %s' %(updown,rPt0)
684 nmohr 1.27 if applyRegression:
685     hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
686     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
687     hJet_MtArray[0][0] = hJ0.Mt()
688     hJet_MtArray[1][0] = hJ1.Mt()
689     hJet_EtArray[0][0] = hJ0.Et()
690     hJet_EtArray[1][0] = hJ1.Et()
691 nmohr 1.37 for key in regVars:
692     var = regDict[key]
693     if key == 'Jet_pt' or key == 'Jet_e' or key == 'hJet_pt' or key == 'hJet_e' or key == 'Jet_ptRaw' or key =='VHbb::evalEt(hJet_pt,hJet_eta,hJet_phi,hJet_e)' or key =='VHbb::evalMt(hJet_pt,hJet_eta,hJet_phi,hJet_e)' or key == 'VHbb::evalJERBias(hJet_ptRaw,hJet_genPt,hJet_eta)':
694     if key == 'Jet_ptRaw':
695     hJet_ptRawArray[0][0] = hJet_ptRaw0*(1+variation*hJet_JECUnc0)
696     hJet_ptRawArray[1][0] = hJet_ptRaw1*(1+variation*hJet_JECUnc1)
697     elif key == 'VHbb::evalJERBias(hJet_ptRaw,hJet_genPt,hJet_eta)':
698     theVars0[key][0] = hJet_ptRaw0*(1+variation*hJet_JECUnc0)
699     theVars1[key][0] = hJet_ptRaw1*(1+variation*hJet_JECUnc1)
700     elif var == 'Jet_pt' or var == 'hJet_pt[0]' or var == 'hJet_pt[1]' :
701     theVars0[key][0] = rPt0
702     theVars1[key][0] = rPt1
703     elif var == 'Jet_e' or var == 'hJet_e':
704     theVars0[key][0] = rE0
705     theVars1[key][0] = rE1
706     elif var == 'hJet_e':
707     theVars0[key][0] = rE0
708     theVars1[key][0] = rE1
709     elif key == 'VHbb::evalEt(hJet_pt,hJet_eta,hJet_phi,hJet_e)':
710     theVars0[key][0] = hJ0.Et()
711     theVars1[key][0] = hJ1.Et()
712     elif key == 'VHbb::evalMt(hJet_pt,hJet_eta,hJet_phi,hJet_e)':
713     theVars0[key][0] = hJ0.Mt()
714     theVars1[key][0] = hJ1.Mt()
715 nmohr 1.27 rPt0 = max(0.0001,readerJet0.EvaluateRegression( "jet0Regression" )[0])
716     rPt1 = max(0.0001,readerJet1.EvaluateRegression( "jet1Regression" )[0])
717 nmohr 1.37 #print 'JES reg: %s' %rPt0
718 nmohr 1.27 rE0 = hJet_e0*rPt0/hJet_pt0
719     rE1 = hJet_e1*rPt1/hJet_pt1
720 nmohr 1.1 hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
721     hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
722     #Fill
723     if updown == 'up':
724     hJet_pt_JES_up[0]=rPt0
725     hJet_pt_JES_up[1]=rPt1
726     hJet_e_JES_up[0]=rE0
727     hJet_e_JES_up[1]=rE1
728     H_JES[0]=(hJ0+hJ1).M()
729     H_JES[2]=(hJ0+hJ1).Pt()
730 nmohr 1.37 HVMass_JES_up[0] = (hJ0+hJ1+vect).M()
731 nmohr 1.1 if updown == 'down':
732     hJet_pt_JES_down[0]=rPt0
733     hJet_pt_JES_down[1]=rPt1
734     hJet_e_JES_down[0]=rE0
735     hJet_e_JES_down[1]=rE1
736     H_JES[1]=(hJ0+hJ1).M()
737     H_JES[3]=(hJ0+hJ1).Pt()
738 nmohr 1.37 HVMass_JES_down[0] = (hJ0+hJ1+vect).M()
739    
740     angleHB_JER_up[0]=fAngleHB_JER_up.EvalInstance()
741     angleHB_JER_down[0]=fAngleHB_JER_down.EvalInstance()
742     angleHB_JES_up[0]=fAngleHB_JES_up.EvalInstance()
743     angleHB_JES_down[0]=fAngleHB_JES_down.EvalInstance()
744     angleZZS[0]=fAngleZZS.EvalInstance()
745     angleZZS_JER_up[0]=fAngleZZS_JER_up.EvalInstance()
746     angleZZS_JER_down[0]=fAngleZZS_JER_down.EvalInstance()
747     angleZZS_JES_up[0]=fAngleZZS_JES_up.EvalInstance()
748     angleZZS_JES_down[0]=fAngleZZS_JES_down.EvalInstance()
749 nmohr 1.1
750     newtree.Fill()
751    
752 nmohr 1.37 print 'Exit loop'
753 nmohr 1.1 newtree.AutoSave()
754 nmohr 1.37 print 'Save'
755 nmohr 1.1 output.Close()
756 nmohr 1.37 print 'Close'