ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C
Revision: 1.41
Committed: Wed Oct 19 10:36:31 2011 UTC (13 years, 6 months ago) by buchmann
Content type: text/plain
Branch: MAIN
CVS Tags: Honeypot, cbaf_2p1ifb
Changes since 1.40: +2 -4 lines
Log Message: 
fixing an issue that appeared for 2 points in mSUGRA space

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4 buchmann 1.31 #include <algorithm>
5     #include <cmath>
6 buchmann 1.1
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 buchmann 1.2 string geqleq;
40     string mcjzbexpression;
41 buchmann 1.4 bool automatized=false;//if we're running this fully automatized we don't want each function to flood the screen
42 buchmann 1.2
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 buchmann 1.1
49 buchmann 1.6 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 buchmann 1.1 //______________________________________________________________________________
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 buchmann 1.7 //____________________________________________________________________________________
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) {
89     iplot++;
90     int count=iplot;
91     // Define new histogram
92     string hname=GetNumericHistoName();
93     TH1F* hossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
94     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSSF,"goff");
95     hname=GetNumericHistoName();
96     TH1F* hossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
97     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSSF,"goff");
98    
99     hname=GetNumericHistoName();
100     TH1F* hosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
101     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSOF,"goff");
102     hname=GetNumericHistoName();
103     TH1F* hosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
104     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
105    
106     hname=GetNumericHistoName();
107     TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
108     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
109     hname=GetNumericHistoName();
110     TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
111     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
112    
113     hname=GetNumericHistoName();
114     TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
115     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
116     hname=GetNumericHistoName();
117     TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
118     events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
119    
120     float obs = hossfp->Integral();
121     float pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
122    
123     delete hossfp,hossfn,hosofp,hosofn;
124     delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
125     return obs-pred;
126     }
127    
128 buchmann 1.1
129     //____________________________________________________________________________________
130     // Efficiency plot
131     TH1F* plotEff(TTree* events, TCut kbase, TString informalname) {
132     iplot++;
133     int count=iplot;
134     // Define new histogram
135     char hname[30]; sprintf(hname,"hJzbEff%d",count);
136     TH1F* hJzbEff = new TH1F(hname,"JZB selection efficiency ; JZB (GeV/c); Efficiency",
137     nBins,jzbMin,jzbMax);
138     Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
139    
140 buchmann 1.6 events->Draw(mcjzbexpression.c_str(),"genJZB>-400"&&kbase,"goff");
141 buchmann 1.1 Float_t maxEff = events->GetSelectedRows();
142 buchmann 1.5 if(verbose>0) dout << hname << " (" << informalname <<") " << maxEff << std::endl;
143 buchmann 1.1
144 buchmann 1.5 if(verbose>0) dout << "JZB max = " << jzbMax << std::endl;
145 buchmann 1.1 // Loop over steps to get efficiency curve
146     char cut[256];
147     for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
148 buchmann 1.6 sprintf(cut,"genJZB>%3f",jzbMin+iBin*step);
149 buchmann 1.2 events->Draw(mcjzbexpression.c_str(),TCut(cut)&&kbase,"goff");
150 buchmann 1.1 Float_t eff = static_cast<Float_t>(events->GetSelectedRows())/maxEff;
151 buchmann 1.5 // dout << "COUCOU " << __LINE__ << std::endl;
152 buchmann 1.1 hJzbEff->SetBinContent(iBin+1,eff);
153     hJzbEff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/maxEff));
154     }
155     return hJzbEff;
156    
157    
158     }
159    
160    
161     //________________________________________________________________________________________
162 pablom 1.23 // Master Formula
163     void master_formula(std::vector<float> eff, float &errHi, float &errLo) {
164    
165     float x0 = eff[0];
166     float deltaPos = 0, deltaNeg = 0;
167     for(int k = 0; k < (eff.size()-1)/2; k++) {
168     float xneg = eff[2*k+2];
169     float xpos = eff[2*k+1];
170     if(xpos-x0>0 || xneg-x0>0) {
171     if(xpos-x0 > xneg-x0) {
172     deltaPos += (xpos-x0)*(xpos-x0);
173     } else {
174     deltaPos += (xneg-x0)*(xneg-x0);
175     }
176     }
177     if(x0-xpos>0 || x0-xneg>0) {
178     if(x0-xpos > x0-xneg) {
179     deltaNeg += (xpos-x0)*(xpos-x0);
180     } else {
181     deltaNeg += (xneg-x0)*(xneg-x0);
182     }
183     }
184     }
185     errHi = sqrt(deltaPos);
186     errLo = sqrt(deltaNeg);
187    
188     }
189    
190    
191     //________________________________________________________________________________________
192     // Get normalization factor for the PDFs
193     float get_norm_pdf_factor(TTree *events, int k) {
194    
195 buchmann 1.25 TH1F *haux = new TH1F("haux", "", 10000, 0, 5);
196 pablom 1.23 char nameVar[20];
197     sprintf(nameVar, "pdfW[%d]", k);
198     events->Project("haux", nameVar);
199 buchmann 1.24 float thisW = haux->Integral();
200     events->Project("haux", "pdfW[0]");
201     float normW = haux->Integral();
202    
203     float factor=thisW/normW;
204 pablom 1.23
205     delete haux;
206    
207     return factor;
208    
209     }
210    
211    
212    
213     //________________________________________________________________________________________
214 buchmann 1.1 // Pile-up efficiency
215 buchmann 1.10 float pileup(TTree *events, bool requireZ, string informalname, string addcut="",Float_t myJzbMax = 140. ) {
216 buchmann 1.1 nBins = 16;
217     jzbMax = myJzbMax;
218    
219     // Acceptance cuts
220 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");
221 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
222    
223 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
224 buchmann 1.1 TH1F* hLM4 = plotEff(events,kbase,informalname);
225     hLM4->SetMinimum(0.);
226    
227     // Nominal function
228     TF1* func = new TF1("func","0.5*TMath::Erfc([0]*x-[1])",jzbMin,jzbMax);
229     func->SetParameter(0,0.03);
230     func->SetParameter(1,0.);
231     hLM4->Fit(func,"Q");
232    
233     // Pimped-up function
234     TF1* funcUp = (TF1*)func->Clone();
235     funcUp->SetParameter( 0., func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
236 buchmann 1.5 if(!automatized) dout << " PU: " << funcUp->Eval(jzbSel) << " " << func->Eval(jzbSel)
237 buchmann 1.1 << "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
238    
239 buchmann 1.18 return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel);
240 buchmann 1.1
241     }
242    
243     //____________________________________________________________________________________
244 buchmann 1.6 // Effect of peak shifting
245 buchmann 1.10 void PeakError(TTree *events,float &result, string mcjzb, float peakerr,string addcut="") {
246 buchmann 1.6 TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
247     TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
248     TString peakcentral("("+TString(mcjzb)+")"+geq_or_leq()+TString(any2string(jzbSel)));
249     TString npeakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
250     TString npeakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
251     TString npeakcentral("("+TString(mcjzb)+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
252     nBins = 1;
253     string informalname="PeakErrorCalculation";
254     float resup,resdown,rescent;
255     for(int i=0;i<3;i++) {
256     string poscut,negcut;
257     if(i==0) {
258     poscut=peakcentral;
259     negcut=npeakcentral;
260     } else if(i==1) {
261     poscut=peakdown;
262     negcut=npeakdown;
263     } else if(i==2) {
264     poscut=peakup;
265     negcut=npeakup;
266     }
267 buchmann 1.10 float res;
268     if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
269     else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
270 buchmann 1.7 if(i==0) rescent=res;
271     else if(i==1) resdown=res;
272     else if(i==2) resup=res;
273 buchmann 1.6 }
274 buchmann 1.41 if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/(float)rescent);
275     else result=(TMath::Abs(rescent-resdown)/(float)rescent);
276 buchmann 1.6 }
277    
278     //____________________________________________________________________________________
279 buchmann 1.1 // Total selection efficiency (MC)
280 buchmann 1.33 //returns the efficiency WITHOUT signal contamination, and the result and resulterr contain the result and the corresponding error
281     Value MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
282 buchmann 1.40 if(!events) {
283     write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
284     result=0;resulterr=0;
285     return Value(0,0);
286     }
287 buchmann 1.1 char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
288     // All acceptance cuts at gen. level
289 buchmann 1.14 //TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
290     TCut kbase("");
291 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
292 buchmann 1.9 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
293 buchmann 1.1 // Corresponding reco. cuts
294 buchmann 1.11 TCut ksel("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
295 buchmann 1.15 TCut ksel2("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
296 pablom 1.23 TCut posSide = kbase&&ksel;
297     TCut negSide = kbase&&ksel2;
298     string sposSide(posSide);
299     string snegSide(negSide);
300     char var[20];
301     sprintf(var, "pdfW[%d]", k);
302     string svar(var);
303     string newPosSide = "(" + sposSide + ")*" + svar;
304     string newNegSide = "(" + snegSide + ")*" + svar;
305    
306 buchmann 1.24 TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
307 buchmann 1.37 if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newPosSide.c_str(),"goff");
308     else events->Draw((mcjzbexpression+">>effh").c_str(), sposSide.c_str(),"goff");
309 buchmann 1.24 Float_t sel = effh->Integral();
310 buchmann 1.21 Float_t nsel=0;
311     if(ConsiderSignalContaminationForLimits) {
312 buchmann 1.38 if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
313     else events->Draw((mcjzbexpression+">>effh").c_str(), snegSide.c_str(),"goff");
314 buchmann 1.24 nsel = effh->Integral();
315 buchmann 1.21 }
316 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.
317 buchmann 1.26 float normFactor = 1;
318     if(k>=0) get_norm_pdf_factor(events, k);
319 pablom 1.23 sel = sel/normFactor;
320     nsel = nsel/normFactor;
321    
322 buchmann 1.13 // events->Draw(mcjzbexpression.c_str(),kbase,"goff");
323     // Float_t tot = events->GetSelectedRows();
324     Float_t tot = Neventsinfile;
325 buchmann 1.1
326 buchmann 1.33 Value result_wo_signalcont;
327    
328 buchmann 1.21 if(ConsiderSignalContaminationForLimits) {
329     result=(sel-nsel)/tot;
330     resulterr=(1.0/tot)*TMath::Sqrt(sel+nsel+(sel-nsel)*(sel-nsel)/tot);
331 buchmann 1.33 result_wo_signalcont=Value(sel/tot,TMath::Sqrt(sel/tot*(1+sel/tot)/tot));
332 buchmann 1.21 } else {//no signal contamination considered:
333     result=(sel)/tot;
334     resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
335 buchmann 1.33 result_wo_signalcont=Value(result,resulterr);
336 buchmann 1.21 }
337 buchmann 1.39 if(!automatized && k>0 ) dout << "PDF assessment (" << k << ") : ";
338 buchmann 1.24 if(!automatized) dout << " MC efficiency: " << result << "+-" << resulterr << " ( JZB>" << jzbSel << " : " << sel << " , JZB<-" << jzbSel << " : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
339     delete effh;
340 buchmann 1.34 return result_wo_signalcont;
341 buchmann 1.1 }
342    
343 buchmann 1.29
344     //____________________________________________________________________________________
345 buchmann 1.30 // Selection efficiency for one process (MC)
346 buchmann 1.29 vector<float> processMCefficiency(TTree *events,string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
347     vector<float> process_efficiencies;
348     for(int iprocess=0;iprocess<=10;iprocess++) {
349     float this_process_efficiency,efferr;
350     stringstream addcutplus;
351     addcutplus<<addcut<<"&&(process=="<<iprocess<<")";
352     MCefficiency(events,this_process_efficiency, efferr,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
353     process_efficiencies.push_back(this_process_efficiency);
354     }
355     return process_efficiencies;
356     }
357    
358    
359 buchmann 1.10 void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, bool requireZ, string addcut="") {
360 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");
361 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
362 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
363 buchmann 1.1 TH1F* hLM4 = plotEff(events,kbase,informalname);
364     Int_t bin = hLM4->FindBin(jzbSel); // To get the error
365 buchmann 1.8 jzbeff=Interpolate(jzbSel,hLM4);
366     jzbefferr=hLM4->GetBinError(bin);
367 buchmann 1.5 if(!automatized) dout << " Efficiency at JZB==" << jzbSel << std::endl;
368 buchmann 1.8 if(!automatized) dout << " " << jzbeff << "+-" << jzbefferr << std::endl;
369 buchmann 1.1 }
370    
371     //________________________________________________________________________
372     // Effect of energy scale on efficiency
373 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 = "" ) {
374 buchmann 1.6 TCut kbase("abs(genMll-91.2)<20&&genZPt>0");
375 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
376 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
377    
378 buchmann 1.1 TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
379     TCut nJets("pfJetGoodNum>2");
380     stringstream down,up;
381     down << "pfJetGoodNum"<<30*(1-syst)<<">=3";
382     up << "pfJetGoodNum"<<30*(1+syst)<<">=3";
383    
384     TCut nJetsP(up.str().c_str());
385     TCut nJetsM(down.str().c_str());
386    
387     if ( !(plotName.Length()>1) ) plotName = informalname;
388    
389     nBins = 1; jzbMin = jzbSel*0.95; jzbMax = jzbSel*1.05;
390     TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname);
391    
392     TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname);
393    
394     TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname);
395    
396     // Dump some information
397     Float_t eff = Interpolate(jzbSel,hist);
398     Float_t effp = Interpolate(jzbSel,histp);
399     Float_t effm = Interpolate(jzbSel,histm);
400 buchmann 1.5 if(!automatized) dout << " Efficiency at JZB==" << jzbSel << std::endl;
401     if(!automatized) dout << " JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
402     if(!automatized) dout << " central: " << eff << std::endl;
403     if(!automatized) dout << " JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
404 buchmann 1.18 jesup=(effp-eff)/eff;
405     jesdown=(effm-eff)/eff;
406 buchmann 1.1 }
407    
408     //________________________________________________________________________
409     // Effect of energy scale on JZB efficiency
410 buchmann 1.10 void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, bool requireZ, string addcut) {
411 buchmann 1.1
412 buchmann 1.6 TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
413 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
414 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
415 buchmann 1.1 TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
416    
417     nBins = 50;
418     jzbMin = 0.5*jzbSel;
419     jzbMax = 2.0*jzbSel;
420    
421     TH1F* hist = plotEff(events,kbase&&ksel,informalname);
422    
423     // Dump some information
424     Float_t eff = Interpolate(jzbSel,hist);
425     Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
426     Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
427 buchmann 1.33 if(!automatized) dout << " efficiency at JZB==" << jzbSel*(1.+systematic) << "(-"<<systematic*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)" << std::endl;
428 buchmann 1.5 if(!automatized) dout << " efficiency at JZB==" << jzbSel << ": " << eff << std::endl;
429 buchmann 1.33 if(!automatized) dout << " efficiency at JZB==" << jzbSel*(1.-systematic) << "(-"<<systematic*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)" << std::endl;
430 buchmann 1.18 up=((effp-eff)/eff);
431     down=((effm-eff)/eff);
432 buchmann 1.1 }
433    
434     //________________________________________________________________________
435     // JZB response (true/reco. vs. true)
436 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 ) {
437 buchmann 1.1
438     jzbMin = 20;
439 buchmann 1.6 TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
440 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
441 buchmann 1.8 if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
442 buchmann 1.1 TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
443    
444 buchmann 1.8 TProfile* hJzbResp = new TProfile("hJzbResp","JZB response ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
445 buchmann 1.1
446 buchmann 1.6 if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZB:genJZB",kbase&&ksel);
447 buchmann 1.1 else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
448    
449     hJzbResp->SetMaximum(1.2);
450     hJzbResp->SetMinimum(0.2);
451     hJzbResp->Fit("pol0","Q");
452     TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
453 buchmann 1.32 if(!fittedfunction) {
454     // in case there are not enough points passing our selection
455     cout << "OOPS response function invalid, assuming 100% error !!!!" << endl;
456     resp=1;
457     resperr=1;
458     } else {
459     resp=fittedfunction->GetParameter(0);
460     resperr=fittedfunction->GetParError(0);
461     if(!automatized) dout << " Response: " << resp << " +/- " << resperr << endl;
462     }
463 buchmann 1.5 delete hJzbResp;
464 buchmann 1.1 }
465    
466    
467 pablom 1.23 //________________________________________________________________________________________
468     // PDF uncertainty
469 buchmann 1.24 float get_pdf_uncertainty(TTree *events, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
470 pablom 1.23 std::vector<float> efficiency;
471 buchmann 1.24 for(int k = 1; k < NPdfs; k++) {
472 pablom 1.23 float result, resulterr;
473     MCefficiency(events, result, resulterr, mcjzb, requireZ, Neventsinfile, addcut, k);
474     efficiency.push_back(result);
475     }
476 buchmann 1.24 float errHi, errLow,err;
477 pablom 1.23 master_formula(efficiency, errHi, errLow);
478 buchmann 1.24 err=errLow;
479     if(errHi>errLow) err=errHi;
480     if(!automatized) dout << " Uncertainty from PDF: " << errLow << " (low) and " << errHi << "(high) ---> Picked " << err << endl;
481     return err;
482 pablom 1.23
483     }
484    
485 buchmann 1.24 int get_npdfs(TTree *events) {
486     int NPDFs;
487     events->SetBranchAddress("NPdfs",&NPDFs);
488     events->GetEntry(1);
489     return NPDFs;
490     }
491    
492 pablom 1.23
493 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) {
494 buchmann 1.1 float JetEnergyScaleUncert=0.1;
495     float JZBScaleUncert=0.1;
496 buchmann 1.2 mcjzbexpression=mcjzb;
497 buchmann 1.28 float triggereff=5.0/100;// in range [0,1]
498 buchmann 1.5 dout << "Trigger efficiency not implemented in this script yet, still using external one" << endl;
499 buchmann 1.18 float leptonseleff=2.0/100;// in range [0,1]
500 buchmann 1.32 leptonseleff=TMath::Sqrt(leptonseleff*leptonseleff+leptonseleff*leptonseleff); // because the 2% is per lepton
501 buchmann 1.5 dout << "Lepton selection efficiency not implemented in this script yet, still using external one" << endl;
502 buchmann 1.1
503 buchmann 1.27 int NPdfs=0;
504     if(ismSUGRA) NPdfs = get_npdfs(events);
505 buchmann 1.29
506 buchmann 1.8 float mceff,mcefferr,jzbeff,jzbefferr;
507 buchmann 1.5 if(!automatized) dout << "MC efficiencies:" << endl;
508 buchmann 1.34 Value mceff_nosigcont = MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut,-1);
509 buchmann 1.35 if(!automatized) cout << " Without signal contamination, we find an efficiency of " << mceff_nosigcont << endl;
510    
511 buchmann 1.32 if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
512 buchmann 1.8 if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
513 buchmann 1.1
514 buchmann 1.39
515     if(!(mceff>0)) write_warning(__FUNCTION__,"Efficiency is zero - the systematics will not be computed!");
516    
517    
518 buchmann 1.6 if(!automatized) dout << "Error from Peak position:" << endl;
519     float sysfrompeak=0;
520 buchmann 1.39 if(mceff>0) PeakError(events,sysfrompeak,mcjzb,peakerror,addcut);
521     else dout << "Not computed." << endl;
522 buchmann 1.6
523 buchmann 1.5 if(!automatized) dout << "Jet energy scale: " << std::endl;
524 buchmann 1.39 float jesup=0,jesdown=0;
525     if(mceff>0) JZBjetScale(events,jesdown,jesup,informalname,requireZ,addcut,JetEnergyScaleUncert);
526     else dout << "Not computed." << endl;
527 buchmann 1.1
528 buchmann 1.5 if(!automatized) dout << "JZB scale: " << std::endl;
529 buchmann 1.39 float scaleup=0,scaledown=0,scalesyst=0;
530     if(mceff>0) doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,requireZ,addcut);
531     else dout << "Not computed." << endl;
532 buchmann 1.1
533 buchmann 1.5 if(!automatized) dout << "JZB response: " << std::endl;
534 buchmann 1.39 float resp=0,resperr=0;
535     if(PlottingSetup::computeJZBresponse && mceff>0) {
536 buchmann 1.32 if(!automatized) dout << "JZB response: " << std::endl;
537     JZBresponse(events,requireZ,resp,resperr,addcut);
538     }
539 buchmann 1.1
540 buchmann 1.5 if(!automatized) dout << "Pileup: " << std::endl;
541 buchmann 1.39 float resolution=0;
542     if(mceff>0) resolution=pileup(events,requireZ,informalname,addcut);
543     else dout << "Not computed." << endl;
544 pablom 1.23
545 buchmann 1.22 float PDFuncert=0;
546 buchmann 1.39 if(!automatized&&mceff>0) dout << "Assessing PDF uncertainty: " << std::endl;
547     if(ismSUGRA&&mceff>0) PDFuncert = get_pdf_uncertainty(events, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
548 pablom 1.23
549 buchmann 1.5 dout << "_______________________________________________" << endl;
550 buchmann 1.10 dout << " SUMMARY FOR " << informalname << " with JZB>" << jzbSel << " (all in %) ";
551     if(addcut!="") dout << "With additional cut: " << addcut;
552     dout << endl;
553 buchmann 1.18 dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl; // in range [0,1]
554     dout << "Trigger efficiency: " << triggereff << endl; // in range [0,1]
555     dout << "Lepton Sel Eff: " << leptonseleff << endl; // in range [0,1]
556     dout << "Jet energy scale: " << jesup << " " << jesdown << endl; // in range [0,1]
557     dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
558     dout << "Resolution : " << resolution << endl; // in range [0,1]
559     dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
560 buchmann 1.24 if(ismSUGRA) dout << "PDF uncertainty : " << PDFuncert << endl; // in range [0,1]
561 buchmann 1.32 if(PlottingSetup::computeJZBefficiency) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
562     if(PlottingSetup::computeJZBresponse)dout << "JZB response : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
563 buchmann 1.1
564 buchmann 1.4 float toterr=0;
565 buchmann 1.18 toterr+=(triggereff)*(triggereff);
566     toterr+=(leptonseleff)*(leptonseleff);
567     if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup*jesup); else toterr+=(jesdown*jesdown);
568     if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup*scaleup); else toterr+=(scaledown*scaledown);
569     toterr+=(resolution*resolution);
570     toterr+=(sysfrompeak*sysfrompeak);
571 buchmann 1.22 if(ismSUGRA) toterr+=(PDFuncert*PDFuncert);
572 buchmann 1.18 dout << "TOTAL SYSTEMATICS: " << TMath::Sqrt(toterr) << " --> " << TMath::Sqrt(toterr)*mceff << endl;
573 buchmann 1.19 float systerr=TMath::Sqrt(toterr)*mceff;
574     toterr=TMath::Sqrt(toterr*mceff*mceff+mcefferr*mcefferr);//also includes stat err!
575    
576     dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*systerr << " (syst) %" << endl;
577     dout << " we thus use the sqrt of the sum of the squares of the stat & syst err, which is : " << 100*toterr << endl;
578 buchmann 1.36 dout << "_______________________________________________" << endl;
579 buchmann 1.19
580     //Do not modify the lines below or mess with the order; this order is expected by all limit calculating functions!
581 buchmann 1.4 vector<float> res;
582     res.push_back(jzbSel);
583     res.push_back(mceff);
584     res.push_back(mcefferr);
585     res.push_back(toterr);
586     res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
587 buchmann 1.18 if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
588     if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
589     res.push_back(fabs(resolution));
590 buchmann 1.35 res.push_back(mceff_nosigcont.getValue());
591     res.push_back(mceff_nosigcont.getError());
592     if(ismSUGRA) res.push_back(PDFuncert);
593 buchmann 1.34 results.push_back(res);
594 buchmann 1.1 }
595    
596 buchmann 1.8 vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
597 buchmann 1.4 automatized=true;
598     vector< vector<float> > systematics;
599 buchmann 1.1 for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
600 buchmann 1.5 dout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
601 buchmann 1.1 for(int ibin=0;ibin<bins.size();ibin++) {
602     jzbSel=bins[ibin];
603 buchmann 1.2 geqleq="geq";
604 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);
605 buchmann 1.1 }//end of bin loop
606     }//end of signal loop
607 buchmann 1.4 return systematics;
608 buchmann 1.1 }