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 |
|
|
//______________________________________________________________________________
|
48 |
|
|
Double_t Interpolate(Double_t x, TH1 *histo)
|
49 |
|
|
{
|
50 |
|
|
// Given a point x, approximates the value via linear interpolation
|
51 |
|
|
// based on the two nearest bin centers
|
52 |
|
|
// Andy Mastbaum 10/21/08
|
53 |
|
|
// in newer ROOT versions but not in the one I have so I had to work around that ...
|
54 |
|
|
|
55 |
|
|
Int_t xbin = histo->FindBin(x);
|
56 |
|
|
Double_t x0,x1,y0,y1;
|
57 |
|
|
|
58 |
|
|
if(x<=histo->GetBinCenter(1)) {
|
59 |
|
|
return histo->GetBinContent(1);
|
60 |
|
|
} else if(x>=histo->GetBinCenter(histo->GetNbinsX())) {
|
61 |
|
|
return histo->GetBinContent(histo->GetNbinsX());
|
62 |
|
|
} else {
|
63 |
|
|
if(x<=histo->GetBinCenter(xbin)) {
|
64 |
|
|
y0 = histo->GetBinContent(xbin-1);
|
65 |
|
|
x0 = histo->GetBinCenter(xbin-1);
|
66 |
|
|
y1 = histo->GetBinContent(xbin);
|
67 |
|
|
x1 = histo->GetBinCenter(xbin);
|
68 |
|
|
} else {
|
69 |
|
|
y0 = histo->GetBinContent(xbin);
|
70 |
|
|
x0 = histo->GetBinCenter(xbin);
|
71 |
|
|
y1 = histo->GetBinContent(xbin+1);
|
72 |
|
|
x1 = histo->GetBinCenter(xbin+1);
|
73 |
|
|
}
|
74 |
|
|
return y0 + (x-x0)*((y1-y0)/(x1-x0));
|
75 |
|
|
}
|
76 |
|
|
}
|
77 |
|
|
|
78 |
|
|
|
79 |
|
|
//____________________________________________________________________________________
|
80 |
|
|
// Efficiency plot
|
81 |
|
|
TH1F* plotEff(TTree* events, TCut kbase, TString informalname) {
|
82 |
|
|
iplot++;
|
83 |
|
|
int count=iplot;
|
84 |
|
|
// Define new histogram
|
85 |
|
|
char hname[30]; sprintf(hname,"hJzbEff%d",count);
|
86 |
|
|
TH1F* hJzbEff = new TH1F(hname,"JZB selection efficiency ; JZB (GeV/c); Efficiency",
|
87 |
|
|
nBins,jzbMin,jzbMax);
|
88 |
|
|
Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
|
89 |
|
|
|
90 |
buchmann |
1.2 |
events->Draw(mcjzbexpression.c_str(),"genJZBSel>-400"&&kbase,"goff");
|
91 |
buchmann |
1.1 |
Float_t maxEff = events->GetSelectedRows();
|
92 |
|
|
if(verbose>0) std::cout << hname << " (" << informalname <<") " << maxEff << std::endl;
|
93 |
|
|
|
94 |
|
|
if(verbose>0) std::cout << "JZB max = " << jzbMax << std::endl;
|
95 |
|
|
// Loop over steps to get efficiency curve
|
96 |
|
|
char cut[256];
|
97 |
|
|
for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
|
98 |
|
|
sprintf(cut,"genJZBSel>%3f",jzbMin+iBin*step);
|
99 |
buchmann |
1.2 |
events->Draw(mcjzbexpression.c_str(),TCut(cut)&&kbase,"goff");
|
100 |
buchmann |
1.1 |
Float_t eff = static_cast<Float_t>(events->GetSelectedRows())/maxEff;
|
101 |
|
|
// std::cout << "COUCOU " << __LINE__ << std::endl;
|
102 |
|
|
hJzbEff->SetBinContent(iBin+1,eff);
|
103 |
|
|
hJzbEff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/maxEff));
|
104 |
|
|
}
|
105 |
|
|
return hJzbEff;
|
106 |
|
|
|
107 |
|
|
|
108 |
|
|
}
|
109 |
|
|
|
110 |
|
|
|
111 |
|
|
//________________________________________________________________________________________
|
112 |
|
|
// Pile-up efficiency
|
113 |
|
|
float pileup(TTree *events, string informalname, Float_t myJzbMax = 140. ) {
|
114 |
|
|
nBins = 16;
|
115 |
|
|
jzbMax = myJzbMax;
|
116 |
|
|
|
117 |
|
|
// Acceptance cuts
|
118 |
buchmann |
1.3 |
TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPtSel>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
|
119 |
buchmann |
1.1 |
|
120 |
|
|
TH1F* hLM4 = plotEff(events,kbase,informalname);
|
121 |
|
|
hLM4->SetMinimum(0.);
|
122 |
|
|
|
123 |
|
|
// Nominal function
|
124 |
|
|
TF1* func = new TF1("func","0.5*TMath::Erfc([0]*x-[1])",jzbMin,jzbMax);
|
125 |
|
|
func->SetParameter(0,0.03);
|
126 |
|
|
func->SetParameter(1,0.);
|
127 |
|
|
hLM4->Fit(func,"Q");
|
128 |
|
|
|
129 |
|
|
// Pimped-up function
|
130 |
|
|
TF1* funcUp = (TF1*)func->Clone();
|
131 |
|
|
funcUp->SetParameter( 0., func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
|
132 |
buchmann |
1.4 |
if(!automatized) std::cout << " PU: " << funcUp->Eval(jzbSel) << " " << func->Eval(jzbSel)
|
133 |
buchmann |
1.1 |
<< "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
|
134 |
|
|
|
135 |
|
|
return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100.;
|
136 |
|
|
|
137 |
|
|
}
|
138 |
|
|
|
139 |
|
|
//____________________________________________________________________________________
|
140 |
|
|
// Total selection efficiency (MC)
|
141 |
buchmann |
1.4 |
void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb) {
|
142 |
buchmann |
1.1 |
|
143 |
|
|
char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
|
144 |
|
|
// All acceptance cuts at gen. level
|
145 |
buchmann |
1.3 |
TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
|
146 |
buchmann |
1.1 |
// Corresponding reco. cuts
|
147 |
buchmann |
1.2 |
TCut ksel("abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
|
148 |
buchmann |
1.1 |
|
149 |
buchmann |
1.2 |
events->Draw(mcjzbexpression.c_str(),kbase&&ksel,"goff");
|
150 |
buchmann |
1.1 |
Float_t sel = events->GetSelectedRows();
|
151 |
buchmann |
1.2 |
events->Draw(mcjzbexpression.c_str(),kbase,"goff");
|
152 |
buchmann |
1.1 |
Float_t tot = events->GetSelectedRows();
|
153 |
|
|
|
154 |
buchmann |
1.4 |
result=sel/tot;
|
155 |
|
|
resulterr=TMath::Sqrt(sel/tot*(1-sel/tot)/tot);
|
156 |
|
|
if(!automatized) std::cout << " MC efficiency: " << result << "+-" << resulterr << std::endl;
|
157 |
buchmann |
1.1 |
}
|
158 |
|
|
|
159 |
|
|
float JZBefficiency(TTree *events, string informalname) {
|
160 |
buchmann |
1.3 |
TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
|
161 |
buchmann |
1.1 |
TH1F* hLM4 = plotEff(events,kbase,informalname);
|
162 |
|
|
Int_t bin = hLM4->FindBin(jzbSel); // To get the error
|
163 |
buchmann |
1.4 |
if(!automatized) std::cout << " Efficiency at JZB==" << jzbSel << std::endl;
|
164 |
|
|
if(!automatized) std::cout << " " << Interpolate(jzbSel,hLM4) << "+-" << hLM4->GetBinError(bin) << std::endl;
|
165 |
buchmann |
1.1 |
return -1;
|
166 |
|
|
}
|
167 |
|
|
|
168 |
|
|
//________________________________________________________________________
|
169 |
|
|
// Effect of energy scale on efficiency
|
170 |
|
|
void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
|
171 |
|
|
TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0");
|
172 |
|
|
TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
|
173 |
|
|
TCut nJets("pfJetGoodNum>2");
|
174 |
|
|
stringstream down,up;
|
175 |
|
|
down << "pfJetGoodNum"<<30*(1-syst)<<">=3";
|
176 |
|
|
up << "pfJetGoodNum"<<30*(1+syst)<<">=3";
|
177 |
|
|
|
178 |
|
|
TCut nJetsP(up.str().c_str());
|
179 |
|
|
TCut nJetsM(down.str().c_str());
|
180 |
|
|
|
181 |
|
|
if ( jzbSelection>0 ) jzbSel = jzbSelection;
|
182 |
|
|
|
183 |
|
|
if ( !(plotName.Length()>1) ) plotName = informalname;
|
184 |
|
|
|
185 |
|
|
nBins = 1; jzbMin = jzbSel*0.95; jzbMax = jzbSel*1.05;
|
186 |
|
|
TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname);
|
187 |
|
|
|
188 |
|
|
TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname);
|
189 |
|
|
|
190 |
|
|
TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname);
|
191 |
|
|
|
192 |
|
|
// Dump some information
|
193 |
|
|
Float_t eff = Interpolate(jzbSel,hist);
|
194 |
|
|
Float_t effp = Interpolate(jzbSel,histp);
|
195 |
|
|
Float_t effm = Interpolate(jzbSel,histm);
|
196 |
buchmann |
1.4 |
if(!automatized) std::cout << " Efficiency at JZB==" << jzbSel << std::endl;
|
197 |
|
|
if(!automatized) std::cout << " JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
|
198 |
|
|
if(!automatized) std::cout << " central: " << eff << std::endl;
|
199 |
|
|
if(!automatized) std::cout << " JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
|
200 |
buchmann |
1.1 |
jesup=(effp-eff)/eff*100.;
|
201 |
|
|
jesdown=(effm-eff)/eff*100.;
|
202 |
|
|
}
|
203 |
|
|
|
204 |
|
|
//________________________________________________________________________
|
205 |
|
|
// Effect of energy scale on JZB efficiency
|
206 |
|
|
void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname) {
|
207 |
|
|
|
208 |
buchmann |
1.3 |
TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0&&genNjets>2");
|
209 |
buchmann |
1.1 |
TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
|
210 |
|
|
|
211 |
|
|
nBins = 50;
|
212 |
|
|
jzbMin = 0.5*jzbSel;
|
213 |
|
|
jzbMax = 2.0*jzbSel;
|
214 |
|
|
|
215 |
|
|
TH1F* hist = plotEff(events,kbase&&ksel,informalname);
|
216 |
|
|
|
217 |
|
|
// Dump some information
|
218 |
|
|
Float_t eff = Interpolate(jzbSel,hist);
|
219 |
|
|
Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
|
220 |
|
|
Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
|
221 |
buchmann |
1.4 |
if(!automatized) std::cout << " efficiency at JZB==" << jzbSel*(1.+systematic) << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)" << std::endl;
|
222 |
|
|
if(!automatized) std::cout << " efficiency at JZB==" << jzbSel << ": " << eff << std::endl;
|
223 |
|
|
if(!automatized) std::cout << " efficiency at JZB==" << jzbSel*(1.-systematic) << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)" << std::endl;
|
224 |
buchmann |
1.1 |
up=((effp-eff)/eff)*100;
|
225 |
|
|
down=((effm-eff)/eff)*100;
|
226 |
|
|
}
|
227 |
|
|
|
228 |
|
|
//________________________________________________________________________
|
229 |
|
|
// JZB response (true/reco. vs. true)
|
230 |
|
|
void JZBresponse(TTree *events, bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
|
231 |
|
|
|
232 |
|
|
jzbMin = 20;
|
233 |
buchmann |
1.3 |
TCut kbase("abs(genMllSel-91.2)<20&&genZPtSel>0&&genNjets>2");
|
234 |
buchmann |
1.1 |
TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
|
235 |
|
|
|
236 |
|
|
TProfile* hJzbResp = new TProfile("hJzbResp","JZB response ; JZB true (GeV/c); JZB reco. / JZB true",
|
237 |
|
|
nPeriods, jzbMin, myJzbMax, "" );
|
238 |
|
|
|
239 |
buchmann |
1.2 |
if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZBSel:genJZBSel",kbase&&ksel);
|
240 |
buchmann |
1.1 |
else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
|
241 |
|
|
|
242 |
|
|
hJzbResp->SetMaximum(1.2);
|
243 |
|
|
hJzbResp->SetMinimum(0.2);
|
244 |
|
|
hJzbResp->Fit("pol0","Q");
|
245 |
|
|
TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
|
246 |
buchmann |
1.4 |
if(!automatized) cout << " Response: " << fittedfunction->GetParameter(0) << " +/- " << fittedfunction->GetParError(0) << endl;
|
247 |
buchmann |
1.1 |
}
|
248 |
|
|
|
249 |
|
|
|
250 |
buchmann |
1.4 |
void do_systematics_for_one_file(TTree *events,string informalname, vector<vector<float> > &results,string mcjzb,string datajzb) {
|
251 |
buchmann |
1.1 |
|
252 |
|
|
float JetEnergyScaleUncert=0.1;
|
253 |
|
|
float JZBScaleUncert=0.1;
|
254 |
buchmann |
1.2 |
mcjzbexpression=mcjzb;
|
255 |
buchmann |
1.1 |
|
256 |
|
|
float triggereff=4;//percent!
|
257 |
|
|
cout << "Trigger efficiency not implemented in this script yet, still using external one" << endl;
|
258 |
|
|
float leptonseleff=2;//percent!
|
259 |
|
|
cout << "Lepton selection efficiency not implemented in this script yet, still using external one" << endl;
|
260 |
|
|
|
261 |
|
|
float mceff,mcefferr;
|
262 |
buchmann |
1.4 |
if(!automatized) cout << "MC efficiencies:" << endl;
|
263 |
buchmann |
1.2 |
MCefficiency(events,mceff,mcefferr,mcjzb);
|
264 |
buchmann |
1.1 |
JZBefficiency(events,informalname);
|
265 |
|
|
|
266 |
buchmann |
1.4 |
if(!automatized) std::cout << "Jet energy scale: " << std::endl;
|
267 |
buchmann |
1.1 |
float jesup,jesdown;
|
268 |
|
|
JZBjetScale(events,jesdown,jesup,informalname,JetEnergyScaleUncert);
|
269 |
|
|
|
270 |
buchmann |
1.4 |
if(!automatized) std::cout << "JZB scale: " << std::endl;
|
271 |
buchmann |
1.1 |
float scaleup,scaledown,scalesyst;
|
272 |
|
|
doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname);
|
273 |
|
|
|
274 |
buchmann |
1.4 |
if(!automatized) std::cout << "JZB response: " << std::endl;
|
275 |
buchmann |
1.1 |
JZBresponse(events);
|
276 |
|
|
|
277 |
buchmann |
1.4 |
if(!automatized) std::cout << "Pileup: " << std::endl;
|
278 |
buchmann |
1.1 |
float resolution=pileup(events,informalname);
|
279 |
|
|
|
280 |
|
|
cout << "_______________________________________________" << endl;
|
281 |
|
|
cout << " SUMMARY FOR " << informalname << " with JZB>" << jzbSel << endl;
|
282 |
buchmann |
1.4 |
cout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl;
|
283 |
buchmann |
1.1 |
cout << "Trigger efficiency: " << triggereff << endl;
|
284 |
|
|
cout << "Lepton Sel Eff: " << leptonseleff << endl;
|
285 |
|
|
cout << "For JZB>" << jzbSel << endl;
|
286 |
|
|
cout << "Jet energy scale: " << jesup << " " << jesdown << " --> suggesting: " << Round(0.5*(fabs(jesup)+fabs(jesdown)),1) << endl;
|
287 |
|
|
cout << "JZB Scale Uncert: " << scaledown << " " << scaleup << " --> suggesting: " << Round(0.5*(fabs(scaledown)+fabs(scaleup)),1) << endl;
|
288 |
|
|
cout << "Resolution : " << resolution << endl;
|
289 |
|
|
|
290 |
|
|
|
291 |
buchmann |
1.4 |
float toterr=0;
|
292 |
|
|
toterr+=(triggereff/100)*(triggereff/100);
|
293 |
|
|
toterr+=(leptonseleff/100)*(leptonseleff/100);
|
294 |
|
|
if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup/100)*(jesup/100); else toterr+=(jesdown/100)*(jesdown/100);
|
295 |
|
|
if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup/100)*(scaleup/100); else toterr+=(scaledown/100)*(scaledown/100);
|
296 |
|
|
toterr+=(resolution/100)*(resolution/100);
|
297 |
|
|
toterr=TMath::Sqrt(toterr);
|
298 |
|
|
cout << "FINAL RESULT : " << mceff << " +/- "<< mcefferr << " (stat) +/- " << 100*toterr << " (syst)" << endl;
|
299 |
|
|
cout << " we thus use the sqrt of the sum of the squares which is : " << 100*TMath::Sqrt(mcefferr*mcefferr+(toterr*toterr)) << endl;
|
300 |
|
|
vector<float> res;
|
301 |
|
|
res.push_back(jzbSel);
|
302 |
|
|
res.push_back(mceff);
|
303 |
|
|
res.push_back(mcefferr);
|
304 |
|
|
res.push_back(toterr);
|
305 |
|
|
res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
|
306 |
|
|
|
307 |
|
|
results.push_back(res);
|
308 |
buchmann |
1.1 |
}
|
309 |
|
|
|
310 |
|
|
vector<vector<float> > compute_systematics(string mcjzb, string datajzb, samplecollection &signalsamples, vector<float> bins) {
|
311 |
buchmann |
1.4 |
automatized=true;
|
312 |
|
|
vector< vector<float> > systematics;
|
313 |
buchmann |
1.1 |
for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
|
314 |
|
|
cout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
|
315 |
|
|
for(int ibin=0;ibin<bins.size();ibin++) {
|
316 |
|
|
jzbSel=bins[ibin];
|
317 |
buchmann |
1.2 |
geqleq="geq";
|
318 |
buchmann |
1.4 |
do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].samplename,systematics,mcjzb,datajzb);
|
319 |
buchmann |
1.1 |
}//end of bin loop
|
320 |
|
|
}//end of signal loop
|
321 |
buchmann |
1.4 |
return systematics;
|
322 |
buchmann |
1.1 |
}
|