ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C
Revision: 1.23
Committed: Wed Aug 31 10:21:17 2011 UTC (13 years, 8 months ago) by pablom
Content type: text/plain
Branch: MAIN
Changes since 1.22: +89 -6 lines
Log Message: 
Added calculation of PDF uncertainty.

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4    
5     #include <TMath.h>
6     #include <TColor.h>
7     #include <TPaveText.h>
8     #include <TRandom.h>
9     #include <TF1.h>
10    
11     #ifndef SampleClassLoaded
12     #include "ActiveSamples.C"
13     #endif
14    
15     #ifndef Verbosity
16     #define Verbosity 0
17     #endif
18    
19     #include <TFile.h>
20     #include <TTree.h>
21     #include <TH1.h>
22     #include <TCut.h>
23     #include <TMath.h>
24     #include <TLine.h>
25     #include <TCanvas.h>
26     #include <TProfile.h>
27     #include <TF1.h>
28    
29    
30    
31     Int_t nBins = 100;
32     Float_t jzbMin = -207;
33     Float_t jzbMax = 243;
34     Float_t jzbSel = 100;
35     int iplot=0;
36     int verbose=0;
37 buchmann 1.2 string geqleq;
38     string mcjzbexpression;
39 buchmann 1.4 bool automatized=false;//if we're running this fully automatized we don't want each function to flood the screen
40 buchmann 1.2
41     TString geq_or_leq() {
42     if(geqleq=="geq") return TString(">=");
43     if(geqleq=="leq") return TString("<=");
44     return TString("GEQ_OR_LEQ_ERROR");
45     }
46 buchmann 1.1
47 buchmann 1.6 TString ngeq_or_leq() {
48     if(geqleq=="geq") return TString("<=");
49     if(geqleq=="leq") return TString(">=");
50     return TString("NGEQ_OR_LEQ_ERROR");
51     }
52    
53 buchmann 1.1 //______________________________________________________________________________
54     Double_t Interpolate(Double_t x, TH1 *histo)
55     {
56     // Given a point x, approximates the value via linear interpolation
57     // based on the two nearest bin centers
58     // Andy Mastbaum 10/21/08
59     // in newer ROOT versions but not in the one I have so I had to work around that ...
60    
61     Int_t xbin = histo->FindBin(x);
62     Double_t x0,x1,y0,y1;
63    
64     if(x<=histo->GetBinCenter(1)) {
65     return histo->GetBinContent(1);
66     } else if(x>=histo->GetBinCenter(histo->GetNbinsX())) {
67     return histo->GetBinContent(histo->GetNbinsX());
68     } else {
69     if(x<=histo->GetBinCenter(xbin)) {
70     y0 = histo->GetBinContent(xbin-1);
71     x0 = histo->GetBinCenter(xbin-1);
72     y1 = histo->GetBinContent(xbin);
73     x1 = histo->GetBinCenter(xbin);
74     } else {
75     y0 = histo->GetBinContent(xbin);
76     x0 = histo->GetBinCenter(xbin);
77     y1 = histo->GetBinContent(xbin+1);
78     x1 = histo->GetBinCenter(xbin+1);
79     }
80     return y0 + (x-x0)*((y1-y0)/(x1-x0));
81     }
82     }
83    
84 buchmann 1.7 //____________________________________________________________________________________
85     // Plotting with all contributions, i.e. sidebands, peak, osof,ossf ... (for a systematic)
86     float allcontributionsplot(TTree* events, TCut kBaseCut, TCut kMassCut, TCut kSidebandCut, TCut JZBPosCut, TCut JZBNegCut) {
87     iplot++;
88     int count=iplot;
89     // Define new histogram
90     string hname=GetNumericHistoName();
91     TH1F* hossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
92     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSSF,"goff");
93     hname=GetNumericHistoName();
94     TH1F* hossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
95     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSSF,"goff");
96    
97     hname=GetNumericHistoName();
98     TH1F* hosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
99     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSOF,"goff");
100     hname=GetNumericHistoName();
101     TH1F* hosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
102     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
103    
104     hname=GetNumericHistoName();
105     TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
106     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
107     hname=GetNumericHistoName();
108     TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
109     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
110    
111     hname=GetNumericHistoName();
112     TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
113     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
114     hname=GetNumericHistoName();
115     TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
116     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
117    
118     float obs = hossfp->Integral();
119     float pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
120    
121     delete hossfp,hossfn,hosofp,hosofn;
122     delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
123     return obs-pred;
124     }
125    
126 buchmann 1.1
127     //____________________________________________________________________________________
128     // Efficiency plot
129     TH1F* plotEff(TTree* events, TCut kbase, TString informalname) {
130     iplot++;
131     int count=iplot;
132     // Define new histogram
133     char hname[30]; sprintf(hname,"hJzbEff%d",count);
134     TH1F* hJzbEff = new TH1F(hname,"JZB selection efficiency ; JZB (GeV/c); Efficiency",
135     nBins,jzbMin,jzbMax);
136     Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
137    
138 buchmann 1.6 events->Draw(mcjzbexpression.c_str(),"genJZB>-400"&&kbase,"goff");
139 buchmann 1.1 Float_t maxEff = events->GetSelectedRows();
140 buchmann 1.5 if(verbose>0) dout << hname << " (" << informalname <<") " << maxEff << std::endl;
141 buchmann 1.1
142 buchmann 1.5 if(verbose>0) dout << "JZB max = " << jzbMax << std::endl;
143 buchmann 1.1 // Loop over steps to get efficiency curve
144     char cut[256];
145     for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
146 buchmann 1.6 sprintf(cut,"genJZB>%3f",jzbMin+iBin*step);
147 buchmann 1.2 events->Draw(mcjzbexpression.c_str(),TCut(cut)&&kbase,"goff");
148 buchmann 1.1 Float_t eff = static_cast<Float_t>(events->GetSelectedRows())/maxEff;
149 buchmann 1.5 // dout << "COUCOU " << __LINE__ << std::endl;
150 buchmann 1.1 hJzbEff->SetBinContent(iBin+1,eff);
151     hJzbEff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/maxEff));
152     }
153     return hJzbEff;
154    
155    
156     }
157    
158    
159     //________________________________________________________________________________________
160 pablom 1.23 // Master Formula
161     void master_formula(std::vector<float> eff, float &errHi, float &errLo) {
162    
163     float x0 = eff[0];
164     float deltaPos = 0, deltaNeg = 0;
165     for(int k = 0; k < (eff.size()-1)/2; k++) {
166     float xneg = eff[2*k+2];
167     float xpos = eff[2*k+1];
168     if(xpos-x0>0 || xneg-x0>0) {
169     if(xpos-x0 > xneg-x0) {
170     deltaPos += (xpos-x0)*(xpos-x0);
171     } else {
172     deltaPos += (xneg-x0)*(xneg-x0);
173     }
174     }
175     if(x0-xpos>0 || x0-xneg>0) {
176     if(x0-xpos > x0-xneg) {
177     deltaNeg += (xpos-x0)*(xpos-x0);
178     } else {
179     deltaNeg += (xneg-x0)*(xneg-x0);
180     }
181     }
182     }
183     errHi = sqrt(deltaPos);
184     errLo = sqrt(deltaNeg);
185    
186     }
187    
188    
189     //________________________________________________________________________________________
190     // Get normalization factor for the PDFs
191     float get_norm_pdf_factor(TTree *events, int k) {
192    
193     TH1F *haux = new TH1F("haux", "", 1000, 0, 15);
194     char nameVar[20];
195     sprintf(nameVar, "pdfW[%d]", k);
196     events->Project("haux", nameVar);
197     float factor = haux->Integral();
198    
199     delete haux;
200    
201     return factor;
202    
203     }
204    
205    
206    
207     //________________________________________________________________________________________
208 buchmann 1.1 // Pile-up efficiency
209 buchmann 1.10 float pileup(TTree *events, bool requireZ, string informalname, string addcut="",Float_t myJzbMax = 140. ) {
210 buchmann 1.1 nBins = 16;
211     jzbMax = myJzbMax;
212    
213     // Acceptance cuts
214 buchmann 1.6 TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
215 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
216    
217 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
218 buchmann 1.1 TH1F* hLM4 = plotEff(events,kbase,informalname);
219     hLM4->SetMinimum(0.);
220    
221     // Nominal function
222     TF1* func = new TF1("func","0.5*TMath::Erfc([0]*x-[1])",jzbMin,jzbMax);
223     func->SetParameter(0,0.03);
224     func->SetParameter(1,0.);
225     hLM4->Fit(func,"Q");
226    
227     // Pimped-up function
228     TF1* funcUp = (TF1*)func->Clone();
229     funcUp->SetParameter( 0., func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
230 buchmann 1.5 if(!automatized) dout << " PU: " << funcUp->Eval(jzbSel) << " " << func->Eval(jzbSel)
231 buchmann 1.1 << "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
232    
233 buchmann 1.18 return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel);
234 buchmann 1.1
235     }
236    
237     //____________________________________________________________________________________
238 buchmann 1.6 // Effect of peak shifting
239 buchmann 1.10 void PeakError(TTree *events,float &result, string mcjzb, float peakerr,string addcut="") {
240 buchmann 1.6 TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
241     TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
242     TString peakcentral("("+TString(mcjzb)+")"+geq_or_leq()+TString(any2string(jzbSel)));
243     TString npeakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
244     TString npeakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
245     TString npeakcentral("("+TString(mcjzb)+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
246    
247     nBins = 1;
248     string informalname="PeakErrorCalculation";
249     float resup,resdown,rescent;
250     for(int i=0;i<3;i++) {
251     string poscut,negcut;
252     if(i==0) {
253     poscut=peakcentral;
254     negcut=npeakcentral;
255     } else if(i==1) {
256     poscut=peakdown;
257     negcut=npeakdown;
258     } else if(i==2) {
259     poscut=peakup;
260     negcut=npeakup;
261     }
262 buchmann 1.10 float res;
263     if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
264     else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
265 buchmann 1.7 if(i==0) rescent=res;
266     else if(i==1) resdown=res;
267     else if(i==2) resup=res;
268 buchmann 1.6 }
269 buchmann 1.18 if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/rescent);
270     else result=(TMath::Abs(rescent-resdown)/rescent);
271 buchmann 1.6 }
272    
273     //____________________________________________________________________________________
274 buchmann 1.1 // Total selection efficiency (MC)
275 pablom 1.23 void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
276 buchmann 1.1
277     char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
278     // All acceptance cuts at gen. level
279 buchmann 1.14 //TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
280     TCut kbase("");
281 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
282 buchmann 1.9 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
283 buchmann 1.1 // Corresponding reco. cuts
284 buchmann 1.11 TCut ksel("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
285 buchmann 1.15 TCut ksel2("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
286 pablom 1.23 TCut posSide = kbase&&ksel;
287     TCut negSide = kbase&&ksel2;
288     string sposSide(posSide);
289     string snegSide(negSide);
290     char var[20];
291     sprintf(var, "pdfW[%d]", k);
292     string svar(var);
293     string newPosSide = "(" + sposSide + ")*" + svar;
294     string newNegSide = "(" + snegSide + ")*" + svar;
295    
296     events->Draw(mcjzbexpression.c_str(), newPosSide.c_str(),"goff");
297 buchmann 1.1 Float_t sel = events->GetSelectedRows();
298 buchmann 1.21 Float_t nsel=0;
299     if(ConsiderSignalContaminationForLimits) {
300 pablom 1.23 events->Draw(mcjzbexpression.c_str(), newNegSide.c_str(),"goff");
301 buchmann 1.21 nsel = events->GetSelectedRows();
302     }
303 pablom 1.23 //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.
304     float normFactor = get_norm_pdf_factor(events, k);
305     sel = sel/normFactor;
306     nsel = nsel/normFactor;
307    
308 buchmann 1.13 // events->Draw(mcjzbexpression.c_str(),kbase,"goff");
309     // Float_t tot = events->GetSelectedRows();
310     Float_t tot = Neventsinfile;
311 buchmann 1.1
312 buchmann 1.21 if(ConsiderSignalContaminationForLimits) {
313     result=(sel-nsel)/tot;
314     resulterr=(1.0/tot)*TMath::Sqrt(sel+nsel+(sel-nsel)*(sel-nsel)/tot);
315     } else {//no signal contamination considered:
316     result=(sel)/tot;
317     resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
318     }
319 buchmann 1.19 if(!automatized) dout << " MC efficiency: " << result << "+-" << resulterr << " ( JZB>" << jzbSel << " : " << sel << " , JZB<-" << jzbSel << " : " << nsel << " and nevents=" << tot << ")" << std::endl;
320 buchmann 1.1 }
321    
322 buchmann 1.10 void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, bool requireZ, string addcut="") {
323 buchmann 1.6 TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
324 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
325 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
326 buchmann 1.1 TH1F* hLM4 = plotEff(events,kbase,informalname);
327     Int_t bin = hLM4->FindBin(jzbSel); // To get the error
328 buchmann 1.8 jzbeff=Interpolate(jzbSel,hLM4);
329     jzbefferr=hLM4->GetBinError(bin);
330 buchmann 1.5 if(!automatized) dout << " Efficiency at JZB==" << jzbSel << std::endl;
331 buchmann 1.8 if(!automatized) dout << " " << jzbeff << "+-" << jzbefferr << std::endl;
332 buchmann 1.1 }
333    
334     //________________________________________________________________________
335     // Effect of energy scale on efficiency
336 buchmann 1.10 void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname,bool requireZ,string addcut="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
337 buchmann 1.6 TCut kbase("abs(genMll-91.2)<20&&genZPt>0");
338 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
339 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
340    
341 buchmann 1.1 TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
342     TCut nJets("pfJetGoodNum>2");
343     stringstream down,up;
344     down << "pfJetGoodNum"<<30*(1-syst)<<">=3";
345     up << "pfJetGoodNum"<<30*(1+syst)<<">=3";
346    
347     TCut nJetsP(up.str().c_str());
348     TCut nJetsM(down.str().c_str());
349    
350     if ( !(plotName.Length()>1) ) plotName = informalname;
351    
352     nBins = 1; jzbMin = jzbSel*0.95; jzbMax = jzbSel*1.05;
353     TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname);
354    
355     TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname);
356    
357     TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname);
358    
359     // Dump some information
360     Float_t eff = Interpolate(jzbSel,hist);
361     Float_t effp = Interpolate(jzbSel,histp);
362     Float_t effm = Interpolate(jzbSel,histm);
363 buchmann 1.5 if(!automatized) dout << " Efficiency at JZB==" << jzbSel << std::endl;
364     if(!automatized) dout << " JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
365     if(!automatized) dout << " central: " << eff << std::endl;
366     if(!automatized) dout << " JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
367 buchmann 1.18 jesup=(effp-eff)/eff;
368     jesdown=(effm-eff)/eff;
369 buchmann 1.1 }
370    
371     //________________________________________________________________________
372     // Effect of energy scale on JZB efficiency
373 buchmann 1.10 void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, bool requireZ, string addcut) {
374 buchmann 1.1
375 buchmann 1.6 TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
376 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
377 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
378 buchmann 1.1 TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
379    
380     nBins = 50;
381     jzbMin = 0.5*jzbSel;
382     jzbMax = 2.0*jzbSel;
383    
384     TH1F* hist = plotEff(events,kbase&&ksel,informalname);
385    
386     // Dump some information
387     Float_t eff = Interpolate(jzbSel,hist);
388     Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
389     Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
390 buchmann 1.5 if(!automatized) dout << " efficiency at JZB==" << jzbSel*(1.+systematic) << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)" << std::endl;
391     if(!automatized) dout << " efficiency at JZB==" << jzbSel << ": " << eff << std::endl;
392     if(!automatized) dout << " efficiency at JZB==" << jzbSel*(1.-systematic) << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)" << std::endl;
393 buchmann 1.18 up=((effp-eff)/eff);
394     down=((effm-eff)/eff);
395 buchmann 1.1 }
396    
397     //________________________________________________________________________
398     // JZB response (true/reco. vs. true)
399 buchmann 1.10 void JZBresponse(TTree *events, bool requireZ, float &resp, float &resperr, string addcut="",bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
400 buchmann 1.1
401     jzbMin = 20;
402 buchmann 1.6 TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
403 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
404 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
405 buchmann 1.1 TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
406    
407 buchmann 1.8 TProfile* hJzbResp = new TProfile("hJzbResp","JZB response ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
408 buchmann 1.1
409 buchmann 1.6 if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZB:genJZB",kbase&&ksel);
410 buchmann 1.1 else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
411    
412     hJzbResp->SetMaximum(1.2);
413     hJzbResp->SetMinimum(0.2);
414     hJzbResp->Fit("pol0","Q");
415     TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
416 buchmann 1.8 resp=fittedfunction->GetParameter(0);
417     resperr=fittedfunction->GetParError(0);
418     if(!automatized) dout << " Response: " << resp << " +/- " << resperr << endl;
419 buchmann 1.5 delete hJzbResp;
420 buchmann 1.1 }
421    
422    
423 pablom 1.23 //________________________________________________________________________________________
424     // PDF uncertainty
425     float get_pdf_uncertainty(TTree *events, string mcjzb, bool requireZ, int Neventsinfile, string addcut="") {
426    
427     int numberOfPDFs = 44;
428     std::vector<float> efficiency;
429     for(int k = 1; k < numberOfPDFs; k++) {
430     float result, resulterr;
431     MCefficiency(events, result, resulterr, mcjzb, requireZ, Neventsinfile, addcut, k);
432     efficiency.push_back(result);
433     }
434     float errHi, errLow;
435     master_formula(efficiency, errHi, errLow);
436     if(errHi>errLow) return errHi;
437     return errLow;
438    
439     }
440    
441    
442    
443 buchmann 1.22 void do_systematics_for_one_file(TTree *events,int Neventsinfile,string informalname, vector<vector<float> > &results,string mcjzb,string datajzb,float peakerror,bool requireZ=false, string addcut="", bool ismSUGRA=false) {
444 buchmann 1.1
445     float JetEnergyScaleUncert=0.1;
446     float JZBScaleUncert=0.1;
447 buchmann 1.2 mcjzbexpression=mcjzb;
448 buchmann 1.18 float triggereff=4.0/100;// in range [0,1]
449 buchmann 1.5 dout << "Trigger efficiency not implemented in this script yet, still using external one" << endl;
450 buchmann 1.18 float leptonseleff=2.0/100;// in range [0,1]
451 buchmann 1.5 dout << "Lepton selection efficiency not implemented in this script yet, still using external one" << endl;
452 buchmann 1.1
453 buchmann 1.8 float mceff,mcefferr,jzbeff,jzbefferr;
454 buchmann 1.5 if(!automatized) dout << "MC efficiencies:" << endl;
455 buchmann 1.13 MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut);
456 buchmann 1.10 JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
457 buchmann 1.8 if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
458 buchmann 1.1
459 buchmann 1.6 if(!automatized) dout << "Error from Peak position:" << endl;
460     float sysfrompeak=0;
461 buchmann 1.10 PeakError(events,sysfrompeak,mcjzb,peakerror,addcut);
462 buchmann 1.6
463 buchmann 1.5 if(!automatized) dout << "Jet energy scale: " << std::endl;
464 buchmann 1.1 float jesup,jesdown;
465 buchmann 1.10 JZBjetScale(events,jesdown,jesup,informalname,requireZ,addcut,JetEnergyScaleUncert);
466 buchmann 1.1
467 buchmann 1.5 if(!automatized) dout << "JZB scale: " << std::endl;
468 buchmann 1.1 float scaleup,scaledown,scalesyst;
469 buchmann 1.10 doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,requireZ,addcut);
470 buchmann 1.1
471 buchmann 1.5 if(!automatized) dout << "JZB response: " << std::endl;
472 buchmann 1.8 float resp,resperr;
473 buchmann 1.10 JZBresponse(events,requireZ,resp,resperr,addcut);
474 buchmann 1.1
475 buchmann 1.5 if(!automatized) dout << "Pileup: " << std::endl;
476 buchmann 1.10 float resolution=pileup(events,requireZ,informalname,addcut);
477 pablom 1.23
478 buchmann 1.22 float PDFuncert=0;
479 pablom 1.23 if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, mcjzb, requireZ, Neventsinfile, addcut);
480    
481 buchmann 1.5 dout << "_______________________________________________" << endl;
482 buchmann 1.10 dout << " SUMMARY FOR " << informalname << " with JZB>" << jzbSel << " (all in %) ";
483     if(addcut!="") dout << "With additional cut: " << addcut;
484     dout << endl;
485 buchmann 1.18 dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl; // in range [0,1]
486     dout << "Trigger efficiency: " << triggereff << endl; // in range [0,1]
487     dout << "Lepton Sel Eff: " << leptonseleff << endl; // in range [0,1]
488     dout << "Jet energy scale: " << jesup << " " << jesdown << endl; // in range [0,1]
489     dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
490     dout << "Resolution : " << resolution << endl; // in range [0,1]
491     dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
492 buchmann 1.22 dout << "PDF uncertainty : " << PDFuncert << " (not yet included below) " << endl; // in range [0,1]
493 buchmann 1.18 dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
494     dout << "JZB response : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
495 buchmann 1.1
496 buchmann 1.4 float toterr=0;
497 buchmann 1.18 toterr+=(triggereff)*(triggereff);
498     toterr+=(leptonseleff)*(leptonseleff);
499     if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup*jesup); else toterr+=(jesdown*jesdown);
500     if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup*scaleup); else toterr+=(scaledown*scaledown);
501     toterr+=(resolution*resolution);
502     toterr+=(sysfrompeak*sysfrompeak);
503 buchmann 1.22 if(ismSUGRA) toterr+=(PDFuncert*PDFuncert);
504 buchmann 1.18 dout << "TOTAL SYSTEMATICS: " << TMath::Sqrt(toterr) << " --> " << TMath::Sqrt(toterr)*mceff << endl;
505 buchmann 1.19 float systerr=TMath::Sqrt(toterr)*mceff;
506     toterr=TMath::Sqrt(toterr*mceff*mceff+mcefferr*mcefferr);//also includes stat err!
507    
508     dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*systerr << " (syst) %" << endl;
509     dout << " we thus use the sqrt of the sum of the squares of the stat & syst err, which is : " << 100*toterr << endl;
510    
511     //Do not modify the lines below or mess with the order; this order is expected by all limit calculating functions!
512 buchmann 1.4 vector<float> res;
513     res.push_back(jzbSel);
514     res.push_back(mceff);
515     res.push_back(mcefferr);
516     res.push_back(toterr);
517     res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
518 buchmann 1.18 if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
519     if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
520     res.push_back(fabs(resolution));
521 buchmann 1.22 if(ismSUGRA) res.push_back(PDFuncert);
522 buchmann 1.4 results.push_back(res);
523 buchmann 1.1 }
524    
525 buchmann 1.8 vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
526 buchmann 1.4 automatized=true;
527     vector< vector<float> > systematics;
528 buchmann 1.1 for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
529 buchmann 1.5 dout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
530 buchmann 1.1 for(int ibin=0;ibin<bins.size();ibin++) {
531     jzbSel=bins[ibin];
532 buchmann 1.2 geqleq="geq";
533 buchmann 1.13 do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,mcjzb,datajzb,mcpeakerror,requireZ);
534 buchmann 1.1 }//end of bin loop
535     }//end of signal loop
536 buchmann 1.4 return systematics;
537 buchmann 1.1 }