ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/write_regression_systematics.py
Revision: 1.39
Committed: Fri Apr 5 13:20:05 2013 UTC (12 years, 1 month ago) by nmohr
Content type: text/x-python
Branch: MAIN
CVS Tags: lhcp_UnblindFix, hcp_Unblind, lhcp_11April, LHCP_PreAppFixAfterFreeze, HEAD
Changes since 1.38: +11 -3 lines
Log Message:
Write to scratch and copy to storage

File Contents

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