ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/interface/TriggerWeight.h
Revision: 1.22
Committed: Sun Sep 30 08:32:57 2012 UTC (12 years, 7 months ago) by arizzi
Content type: text/plain
Branch: MAIN
CVS Tags: EDMV42_Step2_V8, EDMV42_Step2_V7, EDMV42_Step2_V6, EDMV42_Step2_V5a, EDMV42_Step2_V5, tauCandV42, hbbsubstructDev_11, hbbsubstructDev_10, hbbsubstructDev_9, hbbsubstructDev_8, hbbsubstructDev_7, hbbsubstructDev_6, hbbsubstructDev_5, hbbsubstructDev_4, hbbsubstructDev_3, hbbsubstructDev_2, hbbsubstructDev_1, hbbsubstructDev, EDMV42_Step2_V4a, EDMV42_Step2_V4, EDMV42_Step2_V3, HEAD
Branch point for: V42TauCandidate, hbbsubstructDevPostHCP
Changes since 1.21: +2 -2 lines
Log Message:
add protection for many jets

File Contents

# Content
1 #ifndef TRIGGERWEIGHT_H
2 #define TRIGGERWEIGHT_H
3
4 #include "FWCore/ParameterSet/interface/ProcessDesc.h"
5 #include "FWCore/PythonParameterSet/interface/PythonProcessDesc.h"
6 //#include "VHbbAnalysis/VHbbDataFormats/interface/TriggerZnunuCurve.h"
7 #include "VHbbAnalysis/VHbbDataFormats/interface/MultiThresholdEfficiency.h"
8 #include <TH1F.h>
9 #include <TF1.h>
10 #include <TFile.h>
11 #include <TTree.h>
12 #include <iostream>
13 class TriggerWeight
14 {
15 public:
16 TriggerWeight(const edm::ParameterSet& ana) : combiner2Thr(2), combiner1Thr(1)
17 {
18 tscaleHLTmu=openFile(ana,"hltMuFileName");
19 tscaleIDmu=openFile(ana,"idMuFileName");
20 tscaleHLTele1=openFile(ana,"hltEle1FileName");
21 tscaleHLTele2=openFile(ana,"hltEle2FileName");
22 tscaleHLTele1Aug=openFile(ana,"hltEle1AugFileName");
23 tscaleHLTele2Aug=openFile(ana,"hltEle2AugFileName");
24 tscaleID80Ele=openFile(ana,"idEle80FileName");
25 tscaleID95Ele=openFile(ana,"idEle95FileName");
26 tscaleHLTeleJet1=openFile(ana,"hltJetEle1FileName");
27 tscaleHLTeleJet2=openFile(ana,"hltJetEle2FileName");
28 tscaleRecoEle=openFile(ana,"recoEleFileName");
29 // tscalePFMHTele=openFile(ana,"hltPFMHTEleFileName");
30 tscaleSingleEleMay=openFile(ana,"hltSingleEleMayFileName");
31 tscaleSingleEleV4=openFile(ana,"hltSingleEleV4FileName");
32 tscaleHLTmuOr30=openFile(ana,"hltMuOr30FileName");
33
34 tscaleSingleEle2012Awp95=openFile(ana,"hltSingleEle2012Awp95");
35 tscaleSingleEle2012Awp80=openFile(ana,"hltSingleEle2012Awp80");
36 tscaleSingleMuon2012A=openFile(ana,"hltSingleMuon2012A");
37 tscaleDoubleEle2012A_leg8=openFile(ana,"hltDoubleEle2012A_leg8");
38 tscaleDoubleEle2012A_leg17=openFile(ana,"hltDoubleEle2012A_leg17");
39 tscaleDoubleMuon2012A_leg8=openFile(ana,"hltDoubleMuon2012A_leg8");
40 tscaleDoubleMuon2012A_leg17=openFile(ana,"hltDoubleMuon2012A_leg17");
41
42 tscaleDoubleMuon2012A_dZ=openFile(ana,"hltDoubleMuon2012A_dZ");
43 tscaleDoubleEle2012A_dZ=openFile(ana,"hltDoubleEle2012A_dZ");
44
45 tscaleMuPlusWCandPt2012A_legMu=openFile(ana,"hltMuPlusWCandPt2012A_legMu");
46 tscaleMuPlusWCandPt2012A_legW=openFile(ana,"hltMuPlusWCandPt2012A_legW");
47
48 tscaleMuID2012A=openFile(ana,"idMu2012A");
49 tscaleEleID2012A=openFile(ana,"idEle2012A");
50 tscaleEleID2012Awp80=openFile(ana,"idEle2012Awp80");
51
52 if(tscaleHLTmu == 0 || tscaleIDmu == 0)
53 {
54 std::cout << "ERROR: cannot load Muon Trigger efficiencies" << std::endl;
55 }
56
57
58
59 }
60
61 static TTree * openFile(const edm::ParameterSet& ana, const char * name)
62 {
63 TFile *hltMuFile = new TFile (ana.getParameter<std::string> (name).c_str(),"read");
64 if(hltMuFile) return (TTree*) hltMuFile->Get("tree");
65 else return 0;
66 }
67
68 static std::pair<float,float> efficiencyFromPtEta(float pt1, float eta1, TTree *t)
69 {
70 // std::cout << "here " << t << " pt1 " << pt1 << " eta1 " << eta1 << std::endl;
71 float s1 = 1.,err=1.;
72 std::pair<float,float> r(s1,err);
73 if(!t) return r;
74 float ptMin,ptMax,etaMin,etaMax,scale,error;
75 int count = 0;
76 t->SetBranchAddress("ptMin",&ptMin);
77 t->SetBranchAddress("ptMax",&ptMax);
78 t->SetBranchAddress("etaMin",&etaMin);
79 t->SetBranchAddress("etaMax",&etaMax);
80 t->SetBranchAddress("scale",&scale);
81 t->SetBranchAddress("error",&error);
82 float lastPtBin = 200;
83 /* for(int jentry = 0; jentry < t->GetEntries(); jentry++)
84 {
85 t->GetEntry(jentry);
86 if(ptMax >= lastPtBin) lastPtBin =ptMax;
87 }*/
88 for(int jentry = 0; jentry < t->GetEntries(); jentry++)
89 {
90 t->GetEntry(jentry);
91 if(ptMax==lastPtBin) ptMax=1e99;
92 if((pt1 > ptMin) && (pt1 < ptMax) && (eta1 > etaMin) && (eta1 < etaMax))
93 {
94 s1 = scale;
95 err=error;
96 count++;
97 }
98 }
99
100 if(count == 0)
101 {
102 return r;
103 }
104
105 r.first=s1;
106 r.second = err;
107 return (r);
108 }
109
110 float scaleMuIsoHLT(float pt1, float eta1)
111 {
112 return efficiencyFromPtEta(pt1,eta1,tscaleHLTmu).first;
113 }
114
115
116
117 float scaleMuID(float pt1, float eta1)
118 {
119 return efficiencyFromPtEta(pt1,eta1,tscaleIDmu).first;
120 }
121
122
123 double scaleDoubleEle17Ele8Aug( std::vector<float> pt, std::vector<float> eta )
124 {
125 std::vector< std::vector<float> > allEleWithEffs;
126 for(unsigned int j=0; j< pt.size() && j < 10 ; j++)
127 {
128 std::vector<float> thisEleEffs;
129 thisEleEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTele1Aug).first);
130 thisEleEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTele2Aug).first);
131
132 allEleWithEffs.push_back(thisEleEffs);
133 }
134
135 return combiner2Thr.weight<Trigger1High2Loose>(allEleWithEffs);
136
137 }
138
139
140 double scaleDoubleEle17Ele8( std::vector<float> pt, std::vector<float> eta )
141 {
142 std::vector< std::vector<float> > allEleWithEffs;
143 for(unsigned int j=0; j< pt.size() && j<10; j++)
144 {
145 std::vector<float> thisEleEffs;
146 thisEleEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTele1).first);
147 thisEleEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTele2).first);
148
149 allEleWithEffs.push_back(thisEleEffs);
150 }
151
152 return combiner2Thr.weight<Trigger1High2Loose>(allEleWithEffs);
153
154 }
155 double muId2012A( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleMuID2012A).first;}
156 double eleId2012A( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleEleID2012A).first;}
157 double eleId2012Awp80( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleEleID2012Awp80).first;}
158
159 double scaleSingleEleMay( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleSingleEleMay).first;}
160 double scaleSingleEleV4( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleSingleEleV4).first; }
161 double scaleID80Ele( float pt, float eta) { return efficiencyFromPtEta(pt,eta,tscaleID80Ele).first; }
162 double scaleID95Ele( float pt, float eta) { return efficiencyFromPtEta(pt,eta,tscaleID95Ele).first; }
163 double scaleRecoEle( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleRecoEle).first; }
164 double scalePFMHTEle( float MetPFPt){
165 double weightPFMHTrigger=0.;
166
167 //FIXME: read from file
168 if(MetPFPt>0. && MetPFPt<5.) weightPFMHTrigger=0.305;
169 if(MetPFPt>5. && MetPFPt<10.) weightPFMHTrigger=0.351;
170 if(MetPFPt>10. && MetPFPt<15.) weightPFMHTrigger=0.461;
171 if(MetPFPt>15. && MetPFPt<20.) weightPFMHTrigger=0.572;
172 if(MetPFPt>20. && MetPFPt<25.) weightPFMHTrigger=0.713;
173 if(MetPFPt>25. && MetPFPt<30.) weightPFMHTrigger=0.844;
174 if(MetPFPt>30. && MetPFPt<35.) weightPFMHTrigger=0.914;
175 if(MetPFPt>35. && MetPFPt<40.) weightPFMHTrigger=0.939;
176 if(MetPFPt>40. && MetPFPt<45.) weightPFMHTrigger=0.981;
177 if(MetPFPt>45. && MetPFPt<50.) weightPFMHTrigger=0.982;
178 if(MetPFPt>50. && MetPFPt<60.) weightPFMHTrigger=0.993;
179 if(MetPFPt>60. && MetPFPt<70.) weightPFMHTrigger=0.995;
180 if(MetPFPt>70. && MetPFPt<100.) weightPFMHTrigger=0.995;
181 if(MetPFPt>100.) weightPFMHTrigger=1.;
182 return weightPFMHTrigger;
183 }
184
185
186 double scaleJet30Jet25( std::vector<float> pt, std::vector<float> eta)
187 {
188
189 std::vector< std::vector<float> > allJetsWithEffs;
190 for(unsigned int j=0; j< pt.size() && j < 10 ; j++)
191 {
192 std::vector<float> thisJetEffs;
193 thisJetEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTeleJet1).first);
194 thisJetEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTeleJet2).first);
195 // std::cout << " jet pt " << pt[j] << " eta " << eta[j] << " eff1 " << thisJetEffs[0] << " eff2 " << thisJetEffs[1] << std::endl;
196 allJetsWithEffs.push_back(thisJetEffs);
197
198 }
199 float res = combiner2Thr.weight<Trigger1High2Loose>(allJetsWithEffs);
200 // std::cout << "Result is " << res << std::endl;
201 return res;
202 // return combiner2Thr.weight<Trigger1High2Loose>(allJetsWithEffs);
203
204 }
205 /*
206 TF1 fpt("f","1-exp(-0.157*(x-19.3))", 0., 9999999.);
207
208 MET80:
209 TF1 fmet80("f","1/ (1 + exp( -0.0709 * (x - 100.7)))", 0., 9999999.);
210
211 MET100:
212 TF1 fmet100("f","1/ (1 + exp( -0.0679 * (x - 128.8)))", 0., 9999999.);
213 */
214
215 //LP curve used for MET
216 double scaleMetHLT( double met){
217 return 1. / (1. + ( exp( 0.059486 * ( 123.27 - met ))));
218 }
219
220 //MET80 component of the factorized JET+MET trigger
221 double scaleMET80(double et)
222 {
223 return 1. / (1. + exp( -0.0709 * (et - 100.7)));
224 }
225
226 //MET100 component
227 double scaleMET100(double et)
228 {
229 return 1. / (1. + exp( -0.0679 * (et - 128.8)));
230 }
231
232 //Single jet20 efficiency for MET+2CJet20
233 double jet20efficiency( double pt)
234 {
235 if(pt < 10 ) return 0;
236 return 1. - exp(-0.157*(pt-19.3));
237 }
238
239 //combined 2 jets efficiency out of N jets, using jet20 efficiency curve
240 double scale2CentralJet( std::vector<float> pt, std::vector<float> eta)
241 {
242
243 std::vector< std::vector<float> > allJetsWithEffs;
244 for(unsigned int j=0; j< pt.size() && j < 10; j++)
245 {
246 if(fabs(eta[j]) < 2.5)
247 {
248 std::vector<float> thisJetEffs;
249 thisJetEffs.push_back(jet20efficiency(pt[j]));
250 allJetsWithEffs.push_back(thisJetEffs);
251 }
252
253 }
254
255 return combiner1Thr.weight<Trigger2SingleThr>(allJetsWithEffs);
256 }
257
258 //New MET 150
259 double scaleMET150(double et)
260 {
261 return 1./ (1. + exp( -0.129226 * (et - 156.699)));
262 }
263
264 float scaleMuOr30IsoHLT(float pt1, float eta1)
265 {
266 return efficiencyFromPtEta(pt1,eta1,tscaleHLTmuOr30).first;
267 }
268
269
270 //For 2012A HLT_DiCentralPFJet30_PFMHT80, valid for pfMET > 100 GeV:
271 double scaleDiJet30MHT80_2012A(double x)
272 {
273 if(x<100) return 0;
274 return (1e0 - exp(-0.04197*(x-75.73))) * 0.9721 ;
275 }
276 //For 2012B HLT_DiCentralJetSumpT100_dPhi05_DiCentralPFJet60_25_PFMET100_HBHENoiseCleaned, valid for pfMET > 100 GeV:
277 double scaleSumpT100MET100_2012B(double x)
278 {
279 if(x<100) return 0;
280 return (1e0 - exp(-0.06704*(x-96.84))) * 0.9199 ;
281 }
282 //For 2012A+B HLT_PFMET150, valid for pfMET > 150 GeV:
283 double scalePFMET150_2012AB(double x)
284 {
285 if(x<150) return 0;
286 return (1e0 - exp(-0.07135*(x-147.4))) * 0.9707;
287 }
288
289
290 //For 2012A HLT_PFMET150 OR HLT_DiCentralPFJet30_PFMHT80, valid for pfMET > 100 GeV:
291 double scalePFMET150orDiJetMET_2012A(double x)
292 {
293 if(x<100) return 0;
294 return (1e0 - exp(-0.0412*(x-75.52))) * 0.9772;
295 }
296
297 //For 2012B HLT_PFMET150 OR HLT_DiCentralJetSumpT100_dPhi05_DiCentralPFJet60_25_PFMET100_HBHENoiseCleaned, valid for pfMET > 100 GeV:
298 double scalePFMET150orDiJetMET_2012B(double x)
299 {
300 if(x<100) return 0;
301 return (1e0 - exp(-0.05482*(x-95.59))) * 0.9702;
302 }
303
304
305 //For 2012C HLT_PFMET150 OR HLT_DiCentralJetSumpT100_dPhi05_DiCentralPFJet60_25_PFMET100_HBHENoiseCleaned, valid for pfMET > 100 GeV:
306 double scalePFMET150orDiJetMET_2012C(double x)
307 {
308 if(x<100) return 0;
309 return (1e0 - exp(-0.05627*(x-95.15))) * 0.9659;
310 }
311
312
313
314
315
316 float doubleEle2012A( float pt1, float eta1, float pt2, float eta2)
317 {
318 // std::cout << "di ele" << std::endl;
319 float eff1_17 = efficiencyFromPtEta(pt1,eta1,tscaleDoubleEle2012A_leg17).first;
320 float eff2_17 = efficiencyFromPtEta(pt2,eta2,tscaleDoubleEle2012A_leg17).first;
321 float eff1_8 = efficiencyFromPtEta(pt1,eta1,tscaleDoubleEle2012A_leg8).first;
322 float eff2_8 = efficiencyFromPtEta(pt2,eta2,tscaleDoubleEle2012A_leg8).first;
323 // std::cout << tscaleDoubleEle2012A_dZ << std::endl;
324 float eff_dz = efficiencyFromPtEta(eta1,eta2,tscaleDoubleEle2012A_dZ).first; // despite the name pt,eta is actually eta1,eta2
325
326 return (eff1_17 * eff2_8 + eff2_17 * eff1_8 - eff1_17*eff2_17)*eff_dz;
327
328 }
329 float doubleMuon2012A( float pt1, float eta1, float pt2, float eta2)
330 {
331 // std::cout << "di mu" << std::endl;
332 float eff1_17 = efficiencyFromPtEta(pt1,eta1,tscaleDoubleMuon2012A_leg17).first;
333 float eff2_17 = efficiencyFromPtEta(pt2,eta2,tscaleDoubleMuon2012A_leg17).first;
334 float eff1_8 = efficiencyFromPtEta(pt1,eta1,tscaleDoubleMuon2012A_leg8).first;
335 float eff2_8 = efficiencyFromPtEta(pt2,eta2,tscaleDoubleMuon2012A_leg8).first;
336 // std::cout << tscaleDoubleMuon2012A_dZ << std::endl;
337 float eff_dz = efficiencyFromPtEta(eta1,eta2,tscaleDoubleMuon2012A_dZ).first; // despite the name pt,eta is actually eta1,eta2
338
339 return (eff1_17 * eff2_8 + eff2_17 * eff1_8 - eff1_17*eff2_17)*eff_dz;
340
341 }
342
343
344
345 float muPlusWCandPt2012A_legW( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleMuPlusWCandPt2012A_legW).first;}
346 float muPlusWCandPt2012A_legMu( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleMuPlusWCandPt2012A_legMu).first;}
347 float singleEle2012Awp80( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleSingleEle2012Awp80).first;}
348 float singleEle2012Awp95( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleSingleEle2012Awp95).first;}
349 float singleMuon2012A( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleSingleMuon2012A).first;}
350
351
352 private:
353 TTree * tscaleHLTele1;
354 TTree * tscaleHLTele2;
355 TTree * tscaleHLTeleJet1;
356 TTree * tscaleHLTeleJet2;
357 TTree * tscaleID80Ele;
358 TTree * tscaleID95Ele;
359 TTree * tscaleRecoEle;
360 TTree * tscaleHLTmuOr30;
361
362 TTree * tscaleSingleEle2012Awp95;
363 TTree * tscaleSingleEle2012Awp80;
364 TTree * tscaleSingleMuon2012A;
365 TTree * tscaleDoubleEle2012A_leg8;
366 TTree * tscaleDoubleEle2012A_leg17;
367 TTree * tscaleDoubleMuon2012A_leg8;
368 TTree * tscaleDoubleMuon2012A_leg17;
369 TTree * tscaleMuPlusWCandPt2012A_legMu;
370 TTree * tscaleMuPlusWCandPt2012A_legW;
371 TTree * tscaleDoubleEle2012A_dZ;
372 TTree * tscaleDoubleMuon2012A_dZ;
373 // TTree * tscalePFMHTele;
374 TTree * tscaleSingleEleMay;
375 TTree * tscaleSingleEleV4;
376
377 TTree * tscaleHLTele1Aug;
378 TTree * tscaleHLTele2Aug;
379
380 TTree * tscaleMuID2012A;
381 TTree * tscaleEleID2012A;
382 TTree * tscaleEleID2012Awp80;
383
384 TTree * tscaleHLTmu;
385 TTree * tscaleIDmu;
386 MultiThresholdEfficiency combiner2Thr;
387 MultiThresholdEfficiency combiner1Thr;
388 };
389
390 #endif