ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C
Revision: 1.53
Committed: Tue Nov 22 13:02:09 2011 UTC (13 years, 5 months ago) by fronga
Content type: text/plain
Branch: MAIN
Changes since 1.52: +12 -11 lines
Log Message:
Added JZB efficiency plot

File Contents

# Content
1 #include <iostream>
2 #include <vector>
3 #include <sys/stat.h>
4 #include <algorithm>
5 #include <cmath>
6
7 #include <TMath.h>
8 #include <TColor.h>
9 #include <TPaveText.h>
10 #include <TRandom.h>
11 #include <TF1.h>
12
13 #ifndef SampleClassLoaded
14 #include "ActiveSamples.C"
15 #endif
16
17 #ifndef Verbosity
18 #define Verbosity 0
19 #endif
20
21 #include <TFile.h>
22 #include <TTree.h>
23 #include <TH1.h>
24 #include <TCut.h>
25 #include <TMath.h>
26 #include <TLine.h>
27 #include <TCanvas.h>
28 #include <TProfile.h>
29 #include <TF1.h>
30
31
32
33 Int_t nBins = 100;
34 Float_t jzbMin = -207;
35 Float_t jzbMax = 243;
36 Float_t jzbSel = 100;
37 int iplot=0;
38 int verbose=0;
39 string geqleq;
40 string mcjzbexpression;
41 bool automatized=false;//if we're running this fully automatized we don't want each function to flood the screen
42
43 TString geq_or_leq() {
44 if(geqleq=="geq") return TString(">=");
45 if(geqleq=="leq") return TString("<=");
46 return TString("GEQ_OR_LEQ_ERROR");
47 }
48
49 TString ngeq_or_leq() {
50 if(geqleq=="geq") return TString("<=");
51 if(geqleq=="leq") return TString(">=");
52 return TString("NGEQ_OR_LEQ_ERROR");
53 }
54
55 //______________________________________________________________________________
56 Double_t Interpolate(Double_t x, TH1 *histo)
57 {
58 // Given a point x, approximates the value via linear interpolation
59 // based on the two nearest bin centers
60 // Andy Mastbaum 10/21/08
61 // in newer ROOT versions but not in the one I have so I had to work around that ...
62
63 Int_t xbin = histo->FindBin(x);
64 Double_t x0,x1,y0,y1;
65
66 if(x<=histo->GetBinCenter(1)) {
67 return histo->GetBinContent(1);
68 } else if(x>=histo->GetBinCenter(histo->GetNbinsX())) {
69 return histo->GetBinContent(histo->GetNbinsX());
70 } else {
71 if(x<=histo->GetBinCenter(xbin)) {
72 y0 = histo->GetBinContent(xbin-1);
73 x0 = histo->GetBinCenter(xbin-1);
74 y1 = histo->GetBinContent(xbin);
75 x1 = histo->GetBinCenter(xbin);
76 } else {
77 y0 = histo->GetBinContent(xbin);
78 x0 = histo->GetBinCenter(xbin);
79 y1 = histo->GetBinContent(xbin+1);
80 x1 = histo->GetBinCenter(xbin+1);
81 }
82 return y0 + (x-x0)*((y1-y0)/(x1-x0));
83 }
84 }
85
86 //____________________________________________________________________________________
87 // Plotting with all contributions, i.e. sidebands, peak, osof,ossf ... (for a systematic)
88 float allcontributionsplot(TTree* events, TCut kBaseCut, TCut kMassCut, TCut kSidebandCut, TCut JZBPosCut, TCut JZBNegCut, int flipped) {
89 iplot++;
90 int count=iplot;
91 string locmcjzbexpression=mcjzbexpression;
92 // Define new histogram
93 string hname=GetNumericHistoName();
94 TH1F* hossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
95 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSSF,"goff");
96 hname=GetNumericHistoName();
97 TH1F* hossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
98 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSSF,"goff");
99
100 hname=GetNumericHistoName();
101 TH1F* hosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
102 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSOF,"goff");
103 hname=GetNumericHistoName();
104 TH1F* hosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
105 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
106
107 float obs=0;
108 float pred=0;
109 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
110 if(PlottingSetup::RestrictToMassPeak) {
111 hname=GetNumericHistoName();
112 TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
113 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
114 hname=GetNumericHistoName();
115 TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
116 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
117
118 hname=GetNumericHistoName();
119 TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
120 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
121 hname=GetNumericHistoName();
122 TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
123 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
124
125 obs = hossfp->Integral();
126 pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
127
128 if(flipped>0) {
129 obs = hossfn->Integral();
130 pred= hossfp->Integral() - (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
131 }
132 delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
133 } else {
134 obs = hossfp->Integral();
135 pred= hossfn->Integral() + (hosofp->Integral() - hosofn->Integral());
136 if(flipped>0) {
137 obs = hossfn->Integral();
138 pred= hossfp->Integral() - (hosofp->Integral() - hosofn->Integral());;
139 }
140 }
141
142 delete hossfp,hossfn,hosofp,hosofn;
143 return obs-pred;
144 }
145
146
147 //____________________________________________________________________________________
148 // Efficiency plot
149 TH1F* plotEff(TTree* events, TCut kbase, TString informalname, int flipped) {
150 iplot++;
151 int count=iplot;
152 iplot++;
153 int count2=iplot;
154 // Define new histogram
155 char hname[30]; sprintf(hname,"hJzbEff%d",count);
156 char hname2[30]; sprintf(hname2,"hJzbEff%d",count2);
157 TH1F* hJzbEff = new TH1F(hname,"JZB selection efficiency ; JZB [GeV]; Efficiency",nBins,jzbMin,jzbMax);
158 TH1F* hJzbEff2= new TH1F(hname2,"JZB selection efficiency ; JZB [GeV]; Efficiency",1,-14000,14000);
159 Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
160
161 if(flipped==0) events->Draw((mcjzbexpression+">>"+(string)hname2).c_str(),("genJZB>-400"&&kbase),"goff");
162 else events->Draw(("(-"+mcjzbexpression+")>>"+(string)hname2).c_str(),("genJZB>-400"&&kbase),"goff");
163 Float_t maxEff = hJzbEff2->Integral();
164 if(verbose>0) dout << hname << " (" << informalname <<") " << maxEff << std::endl;
165
166 if(verbose>0) dout << "JZB max = " << jzbMax << std::endl;
167 // Loop over steps to get efficiency curve
168 char cut[256];
169 for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
170 sprintf(cut,"genJZB>%3f",jzbMin+iBin*step);
171 if(flipped==0) events->Draw((mcjzbexpression+">>"+(string)hname2).c_str(),(TCut(cut)&&kbase),"goff");
172 if(flipped>0) events->Draw(("(-"+mcjzbexpression+")>>"+(string)hname2).c_str(),(TCut(cut)&&kbase),"goff");
173 Float_t eff = static_cast<Float_t>(hJzbEff2->Integral())/maxEff;
174 hJzbEff->SetBinContent(iBin+1,eff);
175 hJzbEff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/maxEff));
176 }
177 delete hJzbEff2;
178 return hJzbEff;
179 }
180
181
182
183 //________________________________________________________________________________________
184 // Master Formula
185 void master_formula(std::vector<float> eff, float &errHi, float &errLo) {
186
187 float x0 = eff[0];
188 float deltaPos = 0, deltaNeg = 0;
189 for(int k = 0; k < (eff.size()-1)/2; k++) {
190 float xneg = eff[2*k+2];
191 float xpos = eff[2*k+1];
192 if(xpos-x0>0 || xneg-x0>0) {
193 if(xpos-x0 > xneg-x0) {
194 deltaPos += (xpos-x0)*(xpos-x0);
195 } else {
196 deltaPos += (xneg-x0)*(xneg-x0);
197 }
198 }
199 if(x0-xpos>0 || x0-xneg>0) {
200 if(x0-xpos > x0-xneg) {
201 deltaNeg += (xpos-x0)*(xpos-x0);
202 } else {
203 deltaNeg += (xneg-x0)*(xneg-x0);
204 }
205 }
206 }
207 errHi = sqrt(deltaPos);
208 errLo = sqrt(deltaNeg);
209
210 }
211
212
213 //________________________________________________________________________________________
214 // Get normalization factor for the PDFs
215 float get_norm_pdf_factor(TTree *events, int k, string addcut) {
216
217 TH1F *haux = new TH1F("haux", "", 10000, 0, 5);
218 char nameVar[20];
219 sprintf(nameVar, "pdfW[%d]", k);
220 events->Project("haux", nameVar, addcut.c_str());
221 float thisW = haux->Integral();
222 events->Project("haux", "pdfW[0]");
223 float normW = haux->Integral();
224
225 float factor=thisW/normW;
226
227 delete haux;
228
229 return factor;
230
231 }
232
233
234
235 //________________________________________________________________________________________
236 // Pile-up efficiency
237 float pileup(TTree *events, bool requireZ, string informalname, int flipped, string addcut="",Float_t myJzbMax = 140. ) {
238 nBins = 16;
239 jzbMax = myJzbMax;
240
241 // Acceptance cuts
242 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
243 TCut kbase(PlottingSetup::genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
244 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
245
246 if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
247 TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
248 hLM4->SetMinimum(0.);
249
250 // Nominal function
251 TF1* func = new TF1("func","0.5*TMath::Erfc([0]*x-[1])",jzbMin,jzbMax);
252 func->SetParameter(0,0.03);
253 func->SetParameter(1,0.);
254 hLM4->Fit(func,"Q");
255
256 // Pimped-up function
257 TF1* funcUp = (TF1*)func->Clone();
258 funcUp->SetParameter( 0, func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
259 if(!automatized) dout << " PU: " << funcUp->Eval(jzbSel) << " " << func->Eval(jzbSel)
260 << "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
261
262 return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel);
263
264 }
265
266 //____________________________________________________________________________________
267 // Effect of peak shifting
268 void PeakError(TTree *events,float &result, string mcjzb, float peakerr,int flipped,string addcut="") {
269 //Note: the cut used here is something like (JZBEXPRESSION+(peakerr)>50) without all the other cuts, to increase statistics (particularly for scans)
270 TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
271 TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
272 TString peakcentral("("+TString(mcjzb)+")"+geq_or_leq()+TString(any2string(jzbSel)));
273 TString npeakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
274 TString npeakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
275 TString npeakcentral("("+TString(mcjzb)+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
276 nBins = 1;
277 string informalname="PeakErrorCalculation";
278 float resup,resdown,rescent;
279 for(int i=0;i<3;i++) {
280 string poscut,negcut;
281 if(i==0) {
282 poscut=peakcentral;
283 negcut=npeakcentral;
284 } else if(i==1) {
285 poscut=peakdown;
286 negcut=npeakdown;
287 } else if(i==2) {
288 poscut=peakup;
289 negcut=npeakup;
290 }
291 float res;
292 if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
293 else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
294 if(i==0) rescent=res;
295 else if(i==1) resdown=res;
296 else if(i==2) resup=res;
297 }
298 if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/(float)rescent);
299 else result=(TMath::Abs(rescent-resdown)/(float)rescent);
300 cout << " " << result << endl;
301 }
302
303
304 void MCPartialefficiency(TTree *events,float &result, float &resulterr,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0, int type = 0) {
305 if(!events) {
306 write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
307 result=0;resulterr=0;
308 return;
309 }
310
311 char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
312 // All acceptance cuts at gen. level
313 //TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
314 TCut kbase("");
315 if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
316 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
317 // Corresponding reco. cuts
318
319 TCut acceptance("genPt2 != 0");
320 TCut massId(cutmass&&cutOSSF);
321 TCut njets(cutnJets);
322 TCut jzbp;
323 TCut jzbn;
324 if(flipped==0) {
325 jzbp=TCut((TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
326 jzbn=TCut((TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
327 } else {
328 jzbp=TCut(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
329 jzbn=TCut(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
330 }
331 float ntotal = events->Draw("pt1", kbase, "goff");
332 TCut theCut;
333 switch(type) {
334 case 1:
335 theCut = kbase+acceptance;
336 break;
337 case 2:
338 theCut = kbase+massId;
339 break;
340 case 3:
341 theCut = kbase+massId+njets;
342 break;
343 case 4:
344 theCut = kbase+massId+njets+jzbn;
345 break;
346 default:
347 theCut = kbase+massId+njets+jzbn;
348 break;
349 }
350
351 string stheCut(theCut);
352 char var[20];
353 sprintf(var, "pdfW[%d]", k);
354
355 string svar(var);
356 string newtheCut;
357 if(k>0) newtheCut = "(" + stheCut + ")*" + svar;
358 else newtheCut = "(" + stheCut + ")"; // for k==0 or even k==-1 we don't need to evaluate PDFs
359
360 TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
361 if(k>=0) events->Draw((mcjzbexpression+">>effh").c_str(), newtheCut.c_str(),"goff");
362 else events->Draw((mcjzbexpression+">>effh").c_str(), theCut,"goff");
363 Float_t sel = effh->Integral();
364 Float_t nsel=0;
365 //Corrections due to normalization in the PDF. This has to be applied as well to the number of events in a file if the definition changes at some point.
366 float normFactor = 1;
367 if(k>=0) get_norm_pdf_factor(events, k, addcut);
368 sel = sel/normFactor;
369
370 result=(sel)/ntotal;
371 resulterr=TMath::Sqrt(sel/ntotal*(1+sel/ntotal)/ntotal);
372
373 delete effh;
374 }
375
376 //____________________________________________________________________________________
377 // Total selection efficiency (MC)
378 //returns the efficiency WITHOUT signal contamination, and the result and resulterr contain the result and the corresponding error
379 Value MCefficiency(TTree *events,float &result, float &resulterr, int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
380 if(!events) {
381 write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
382 result=0;
383 resulterr=0;
384 return Value(0,0);
385 }
386
387 char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
388 // All acceptance cuts at gen. level
389 //TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
390 TCut kbase("");
391
392 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
393 if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
394 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
395 // Corresponding reco. cuts
396
397 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
398 TCut ksel;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
399 TCut ksel2;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
400 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
401 if(PlottingSetup::RestrictToMassPeak||!ConsiderSignalContaminationForLimits) {
402 ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
403 ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
404 if(flipped>0) {
405 ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
406 ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
407 }
408 } else {
409 //for off peak analysis we don't use the OSSF condition here yet so we can recycle these two cuts for the em condition!
410 ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
411 ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
412 if(flipped>0) {
413 ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
414 ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
415 }
416 }
417
418 TCut posSide = kbase&&ksel;
419 TCut negSide = kbase&&ksel2;
420 string sposSide(posSide);
421 string snegSide(negSide);
422 char var[20];
423 sprintf(var, "pdfW[%d]", k);
424 if(k==-1) sprintf(var,"1.0");//case in which we don't want to evaluate PDFs
425 string svar(var);
426 string newPosSide = "((id1==id2)&&(" + sposSide + "))*" + svar;
427 string newNegSide = "((id1==id2)&&(" + snegSide + "))*" + svar;
428 string emnewPosSide = "((id1!=id2)&&(" + sposSide + "))*" + svar; // only used for off peak analysis
429 string emnewNegSide = "((id1!=id2)&&(" + snegSide + "))*" + svar; // only used for off peak analysis
430
431 TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
432 if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), TCut(newPosSide.c_str())*PlottingSetup::Weight,"goff");
433 else events->Draw((mcjzbexpression+">>effh").c_str(), TCut((sposSide+"&&(id1==id2)").c_str())*PlottingSetup::Weight,"goff");//the OSSF condition is added for the offpeak analysis, in onpeak case it's there already but doesn't change anything.
434
435 Float_t sel = effh->Integral();
436 Float_t nsel=0;
437
438 ///----------------------------------------------- THIS PART REQUIRES STUDYING! -------------------------
439
440 if(ConsiderSignalContaminationForLimits) {
441 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
442 if(PlottingSetup::RestrictToMassPeak) {
443 events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(newNegSide.c_str()))*PlottingSetup::Weight,"goff");
444 nsel += effh->Integral();
445 } else {
446 events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(newNegSide.c_str()))*PlottingSetup::Weight,"goff");
447 nsel += effh->Integral();
448 events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(emnewPosSide.c_str()))*PlottingSetup::Weight,"goff");
449 nsel += effh->Integral();
450 events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(emnewNegSide.c_str()))*PlottingSetup::Weight,"goff");
451 nsel -= effh->Integral();
452 }
453 }
454
455 //Corrections due to normalization in the PDF. This has to be applied as well to the number of events in a file if the definition changes at some point.
456 float normFactor = 1;
457 if(k>=0) get_norm_pdf_factor(events, k, addcut);
458 sel = sel/normFactor;
459 nsel = nsel/normFactor;
460
461 Float_t tot = Neventsinfile;
462
463 Value result_wo_signalcont;
464
465 if(ConsiderSignalContaminationForLimits) {
466 result=(sel-nsel)/tot;
467 resulterr=(1.0/tot)*TMath::Sqrt(sel+nsel+(sel-nsel)*(sel-nsel)/tot);
468 result_wo_signalcont=Value(sel/tot,TMath::Sqrt(sel/tot*(1+sel/tot)/tot));
469 } else {//no signal contamination considered:
470 result=(sel)/tot;
471 resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
472 result_wo_signalcont=Value(result,resulterr);
473 }
474 if(!automatized && k>0 ) dout << "PDF assessment [" << k << "] : ";
475 if(!automatized) dout << " MC efficiency: " << result << "+-" << resulterr << " ( JZB>" << jzbSel << " : " << sel << " , signal contamination : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
476 delete effh;
477 return result_wo_signalcont;
478 }
479
480
481
482 //____________________________________________________________________________________
483 // Selection efficiency for one process (MC)
484 // not in use anymore.
485 /*
486 vector<float> processMCefficiency(TTree *events,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
487 vector<float> process_efficiencies;
488 for(int iprocess=0;iprocess<=10;iprocess++) {
489 float this_process_efficiency,efferr;
490 stringstream addcutplus;
491 addcutplus<<addcut<<"&&(process=="<<iprocess<<")";
492 MCefficiency(events,this_process_efficiency, efferr,flipped,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
493 process_efficiencies.push_back(this_process_efficiency);
494 }
495 return process_efficiencies;
496 }
497 */
498
499 void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, int flipped, bool requireZ, string addcut="") {
500 TCut kbase(genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
501 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
502 if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
503 TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
504 Int_t bin = hLM4->FindBin(jzbSel); // To get the error
505 jzbeff=Interpolate(jzbSel,hLM4);
506 jzbefferr=hLM4->GetBinError(bin);
507 if(!automatized) dout << " Efficiency at JZB==" << jzbSel << std::endl;
508 if(!automatized) dout << " " << jzbeff << "+-" << jzbefferr << std::endl;
509 }
510
511 //________________________________________________________________________
512 // Effect of energy scale on efficiency
513 void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname, int flipped, bool requireZ,string addcut="",Float_t jzbSelection=-1, TString plotName = "" ) {
514 TCut kbase(genMassCut&&"genZPt>0");
515 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
516 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
517 if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
518
519 TCut ksel(cutmass&&cutOSSF);
520 TCut nJets("pfJetGoodNum>2");
521 stringstream down,up;
522 down << "pfJetGoodNumn1sigma>=3";
523 up << "pfJetGoodNump1sigma>=3";
524
525 TCut nJetsP(up.str().c_str());
526 TCut nJetsM(down.str().c_str());
527
528 if ( !(plotName.Length()>1) ) plotName = informalname;
529
530 nBins = 1; jzbMin = jzbSel*0.95; jzbMax = jzbSel*1.05;
531 TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname,flipped);
532 TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname,flipped);
533 TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname,flipped);
534
535 // Dump some information
536 Float_t eff = Interpolate(jzbSel,hist);
537 Float_t effp = Interpolate(jzbSel,histp);
538 Float_t effm = Interpolate(jzbSel,histm);
539 if(!automatized) dout << " Efficiency at JZB==" << jzbSel << std::endl;
540 if(!automatized) dout << " JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
541 if(!automatized) dout << " central: " << eff << std::endl;
542 if(!automatized) dout << " JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
543 jesup=(effp-eff)/eff;
544 jesdown=(effm-eff)/eff;
545 }
546
547 //________________________________________________________________________
548 // Effect of energy scale on JZB efficiency
549 void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, int flipped, bool requireZ, string addcut) {
550
551 TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
552 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
553 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
554 if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
555 TCut ksel(cutmass&&cutOSSF);
556
557 nBins = 50;
558 jzbMin = 0.5*jzbSel;
559 jzbMax = 2.0*jzbSel;
560
561 TH1F* hist = plotEff(events,kbase&&ksel,informalname,flipped);
562
563 // Dump some information
564 Float_t eff = Interpolate(jzbSel,hist);
565 Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
566 Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
567 if(!automatized) dout << " efficiency at JZB==" << jzbSel*(1.+systematic) << "(-"<<systematic*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)" << std::endl;
568 if(!automatized) dout << " efficiency at JZB==" << jzbSel << ": " << eff << std::endl;
569 if(!automatized) dout << " efficiency at JZB==" << jzbSel*(1.-systematic) << "(-"<<systematic*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)" << std::endl;
570 up=((effp-eff)/eff);
571 down=((effm-eff)/eff);
572 }
573
574 //________________________________________________________________________
575 // JZB response (true/reco. vs. true)
576 void JZBresponse(TTree *events, string name, bool requireZ, float &resp, float &resperr, int flipped, string addcut="", bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
577
578 jzbMin = 20;
579 flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
580 TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
581 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
582 flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
583 if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
584 flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
585 TCut ksel(cutmass&&cutOSSF);
586
587 TProfile* hJzbResp = new TProfile("hJzbResp","JZB response ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
588
589 string locmcjzbexpression=mcjzbexpression;
590 if(flipped>0) locmcjzbexpression="-"+locmcjzbexpression;
591 string possibleminus="";
592 if(flipped>0) possibleminus="-";
593 if (!isMET) events->Project("hJzbResp","("+TString(locmcjzbexpression)+")/("+possibleminus+"genJZB):("+possibleminus+"genJZB)",kbase&&ksel);
594 else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
595
596 hJzbResp->SetMaximum(1.2);
597 hJzbResp->SetMinimum(0.2);
598 hJzbResp->Fit("pol0","Q");
599 TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
600 if(!fittedfunction) {
601 // in case there are not enough points passing our selection
602 cout << "OOPS response function invalid, assuming 100% error !!!!" << endl;
603 resp=1;
604 resperr=1;
605 } else {
606 resp=fittedfunction->GetParameter(0);
607 resperr=fittedfunction->GetParError(0);
608 if(!automatized) dout << " Response: " << resp << " +/- " << resperr << endl;
609 }
610 delete hJzbResp;
611 }
612
613
614 //________________________________________________________________________________________
615 // PDF uncertainty
616 float get_pdf_uncertainty(TTree *events, int flipped, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
617 std::vector<float> efficiency;
618 for(int k = 1; k < NPdfs; k++) {
619 float result, resulterr;
620 Value flipval;
621 MCefficiency(events, result, resulterr, flipped, mcjzb, requireZ, Neventsinfile, addcut, k);
622 efficiency.push_back(result);
623 }
624 float errHi, errLow,err;
625 master_formula(efficiency, errHi, errLow);
626 err=errLow;
627 if(errHi>errLow) err=errHi;
628 if(!automatized) dout << " Uncertainty from PDF: " << errLow << " (low) and " << errHi << "(high) ---> Picked " << err << endl;
629 return err;
630
631 }
632
633 int get_npdfs(TTree *events) {
634 int NPDFs;
635 events->SetBranchAddress("NPdfs",&NPDFs);
636 events->GetEntry(1);
637 return NPDFs;
638 }
639
640
641 void do_systematics_for_one_file(TTree *events,int Neventsinfile,string informalname, vector<vector<float> > &results,int flipped, string mcjzb,string datajzb,float peakerror,bool requireZ=false, string addcut="", bool ismSUGRA=false) {
642 float JZBScaleUncert=0.05;
643 mcjzbexpression=mcjzb;
644 float triggereff=2.0/100;// in range [0,1]
645 dout << "Trigger efficiency not implemented in this script yet, still using external one" << endl;
646 float leptonseleff=2.0/100;// in range [0,1]
647 leptonseleff=TMath::Sqrt(leptonseleff*leptonseleff+leptonseleff*leptonseleff); // because the 2% is per lepton
648 dout << "Lepton selection efficiency not implemented in this script yet, still using external one" << endl;
649
650 int NPdfs=0;
651 if(ismSUGRA) NPdfs = get_npdfs(events);
652
653 float mceff,mcefferr,jzbeff,jzbefferr;
654 if(!automatized) dout << "MC efficiencies:" << endl;
655 Value flipefficiency;
656 Value mceff_nosigcont = MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,addcut,-1);
657 if(!automatized) cout << " Without signal contamination, we find an efficiency of " << mceff_nosigcont << endl;
658
659 if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,flipped,requireZ,addcut);
660 if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
661
662 if(!automatized) dout << "Error from Peak position:";
663 float sysfrompeak=0;
664 PeakError(events,sysfrompeak,mcjzb,peakerror,flipped,addcut);
665
666 if(!automatized) dout << "Jet energy scale (JES): " << std::endl;
667 float jesup,jesdown;
668 JZBjetScale(events,jesdown,jesup,informalname,flipped,requireZ,addcut);
669
670 if(!automatized) dout << "JZB scale: " << std::endl;
671 float scaleup,scaledown,scalesyst;
672 doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
673
674 if(!automatized) dout << "JZB response: " << std::endl;
675 float resp,resperr;
676 if(PlottingSetup::computeJZBresponse) {
677 if(!automatized) dout << "JZB response: " << std::endl;
678 if(!ismSUGRA) JZBresponse(events,informalname,requireZ,resp,resperr,flipped,addcut);
679 }
680
681 if(!automatized) dout << "Pileup: " << std::endl;
682 // float resolution;
683 //resolution=pileup(events,requireZ,informalname,flipped,addcut);
684
685 float PDFuncert=0;
686 if(!automatized) dout << "Assessing PDF uncertainty: " << std::endl;
687 if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
688
689 dout << "_______________________________________________" << endl;
690 dout << " SUMMARY FOR " << informalname << " with JZB>" << jzbSel << " (all in %) ";
691 if(addcut!="") dout << "With additional cut: " << addcut;
692 dout << endl;
693 dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl; // in range [0,1]
694 dout << "Trigger efficiency: " << triggereff << endl; // in range [0,1]
695 dout << "Lepton Sel Eff: " << leptonseleff << endl; // in range [0,1]
696 dout << "Jet energy scale: " << jesup << " " << jesdown << endl; // in range [0,1]
697 dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
698 // dout << "Resolution : " << resolution << endl; // in range [0,1]
699 dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
700 if(ismSUGRA) dout << "PDF uncertainty : " << PDFuncert << endl; // in range [0,1]
701 if(PlottingSetup::computeJZBefficiency) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
702 if(PlottingSetup::computeJZBresponse)dout << "JZB response : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
703
704 float toterr=0;
705 toterr+=(triggereff)*(triggereff);
706 toterr+=(leptonseleff)*(leptonseleff);
707 if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup*jesup); else toterr+=(jesdown*jesdown);
708 if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup*scaleup); else toterr+=(scaledown*scaledown);
709 // toterr+=(resolution*resolution);
710 toterr+=(sysfrompeak*sysfrompeak);
711 if(ismSUGRA) toterr+=(PDFuncert*PDFuncert);
712 dout << "TOTAL SYSTEMATICS: " << TMath::Sqrt(toterr) << " --> " << TMath::Sqrt(toterr)*mceff << endl;
713 float systerr=TMath::Sqrt(toterr)*mceff;
714 toterr=TMath::Sqrt(toterr*mceff*mceff+mcefferr*mcefferr);//also includes stat err!
715
716 dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*systerr << " (syst) %" << endl;
717 dout << " we thus use the sqrt of the sum of the squares of the stat & syst err, which is : " << 100*toterr << endl;
718 dout << "_______________________________________________" << endl;
719
720 //Do not modify the lines below or mess with the order; this order is expected by all limit calculating functions!
721 vector<float> res;
722 res.push_back(jzbSel);
723 res.push_back(mceff);
724 res.push_back(mcefferr);
725 res.push_back(toterr);
726 res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
727 if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
728 if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
729 // res.push_back(fabs(resolution));
730 res.push_back(0.0);
731 res.push_back(mceff_nosigcont.getValue());
732 res.push_back(mceff_nosigcont.getError());
733 if(ismSUGRA) res.push_back(PDFuncert);
734 results.push_back(res);
735 }
736
737 vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, int flipped, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
738 automatized=true;
739 vector< vector<float> > systematics;
740 for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
741 dout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
742 for(int ibin=0;ibin<bins.size();ibin++) {
743 jzbSel=bins[ibin];
744 geqleq="geq";
745 do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,flipped,mcjzb,datajzb,mcpeakerror,requireZ);
746 }//end of bin loop
747 }//end of signal loop
748 return systematics;
749 }