| 1 |
|
#define qcd_cxx |
| 2 |
|
#include "qcd.h" |
| 3 |
|
#include <TH2.h> |
| 4 |
+ |
#include <THStack.h> |
| 5 |
|
#include <TStyle.h> |
| 6 |
|
#include <TCanvas.h> |
| 7 |
|
#include <TMath.h> |
| 14 |
|
#include <TLegend.h> |
| 15 |
|
#include <TPaveStats.h> |
| 16 |
|
|
| 17 |
< |
int nbin = 50; // 200 for 1 and 2 jet bin |
| 17 |
> |
// Histo range |
| 18 |
|
double xMin = 0.0; |
| 19 |
|
double xMax = 2.0; |
| 19 |
– |
|
| 20 |
– |
// Select fitting function: |
| 21 |
– |
TString fitfunc = "landau"; // landau, polX, gaus ... |
| 22 |
– |
// Select New or Old reliso: |
| 23 |
– |
TString isoname = "New"; |
| 20 |
|
// Jet bin: |
| 21 |
|
Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5 |
| 22 |
+ |
double binwidth = 0.1; // 0.02,0.02,0.05,0.1 |
| 23 |
+ |
int nbin = (xMax-xMin)/binwidth; |
| 24 |
+ |
// fit region weight |
| 25 |
+ |
double wxmin = 1.9; |
| 26 |
+ |
double wxmax = 0.7; |
| 27 |
+ |
// Specify fitting function: |
| 28 |
+ |
TString fitfunc = "landau"; // this version is cut for Landau |
| 29 |
+ |
// Select New or Old reliso: |
| 30 |
+ |
TString isoname = "New"; |
| 31 |
|
double d0sigCut = 3.; |
| 32 |
|
|
| 33 |
|
void qcd::Loop() |
| 34 |
|
{ |
| 35 |
|
// In a ROOT session, you can do: |
| 36 |
+ |
// Modify jbin, binwidth, wxmin, wxmax |
| 37 |
|
// Root > .L qcd.C++ |
| 38 |
|
// Root > qcd t |
| 39 |
|
// Root > t.GetEntry(12); // Fill t data members with entry number 12 |
| 66 |
|
// new |
| 67 |
|
AX0 = 0.0; |
| 68 |
|
AXf = 0.053; |
| 69 |
< |
BX0 = 0.15; |
| 70 |
< |
BXf = 0.6; |
| 69 |
> |
BX0 = 0.4; |
| 70 |
> |
BXf = 2.0; |
| 71 |
|
} |
| 66 |
– |
// // old weights |
| 67 |
– |
// double wttbar = 0.0094; |
| 68 |
– |
// double wqcd = 0.3907; |
| 69 |
– |
// double wwjets = 0.0177; |
| 70 |
– |
// double wzjets = 0.041; |
| 72 |
|
|
| 73 |
|
// weights |
| 74 |
|
double wttbar = 0.0101; |
| 77 |
|
double wzjets = 0.078; |
| 78 |
|
|
| 79 |
|
// Na = num of qcd in signal region |
| 80 |
< |
double Na,Nb,Nfac,Ns, Ns2, NsErrP, NsErrM; |
| 81 |
< |
Na = Nb = Nfac = Ns = Ns2 = NsErrP = NsErrM = 0; |
| 80 |
> |
double Na,Nb,Nfac; |
| 81 |
> |
Na = Nb = Nfac = 0; |
| 82 |
|
|
| 83 |
|
double Nttbar, NWjets, NZjets, Nqcd; |
| 84 |
|
Nttbar = NWjets = NZjets = Nqcd = 0; |
| 158 |
|
|
| 159 |
|
gStyle->SetOptFit(01111); |
| 160 |
|
|
| 161 |
< |
TCanvas *cv1 = new TCanvas("cv1","cv1",800,600); |
| 161 |
> |
TCanvas *cv1 = new TCanvas("cv1","cv1",800,600); |
| 162 |
> |
THStack *hs = new THStack("hs","RelIso (stacked)"); |
| 163 |
|
|
| 164 |
+ |
// Take care of scaling |
| 165 |
|
h_RelIso_qcd->Scale(wqcd); |
| 166 |
|
h_RelIso_wjets->Scale(wwjets); |
| 167 |
|
h_RelIso_zjets->Scale(wzjets); |
| 168 |
|
h_RelIso_ttbar->Scale(wttbar); |
| 169 |
< |
|
| 170 |
< |
h_RelIso_ttbar->SetXTitle("Combined Relative Isolation"); |
| 169 |
> |
|
| 170 |
> |
h_RelIso_qcd->SetXTitle("Combined Relative Isolation"); |
| 171 |
> |
h_RelIso_qcd->SetLineColor(40); |
| 172 |
> |
h_RelIso_qcd->SetFillStyle(3018); |
| 173 |
> |
h_RelIso_qcd->SetFillColor(38); |
| 174 |
> |
hs->Add(h_RelIso_qcd); |
| 175 |
> |
h_RelIso_wjets->SetMarkerColor(4); |
| 176 |
> |
h_RelIso_wjets->SetLineColor(4); |
| 177 |
> |
h_RelIso_wjets->SetFillColor(4); |
| 178 |
> |
hs->Add(h_RelIso_wjets); |
| 179 |
> |
h_RelIso_zjets->SetMarkerColor(5); |
| 180 |
> |
h_RelIso_zjets->SetLineColor(5); |
| 181 |
> |
h_RelIso_zjets->SetFillColor(5); |
| 182 |
> |
hs->Add(h_RelIso_zjets); |
| 183 |
|
h_RelIso_ttbar->SetMarkerColor(2); |
| 184 |
< |
h_RelIso_ttbar->SetLineColor(2); |
| 185 |
< |
h_RelIso_ttbar->Draw(); |
| 184 |
> |
h_RelIso_ttbar->SetLineColor(2); |
| 185 |
> |
h_RelIso_ttbar->SetFillColor(2); |
| 186 |
> |
hs->Add(h_RelIso_ttbar); |
| 187 |
> |
hs->Draw(); |
| 188 |
> |
|
| 189 |
> |
// Label histo |
| 190 |
> |
TLegend * leg0 = new TLegend(0.3,0.3,0.6,0.55); |
| 191 |
> |
leg0->AddEntry(h_RelIso_ttbar,"t#bar{t}"); |
| 192 |
> |
leg0->AddEntry(h_RelIso_qcd,"QCD (MuPt15) "); |
| 193 |
> |
leg0->AddEntry(h_RelIso_wjets,"W+jets"); |
| 194 |
> |
leg0->AddEntry(h_RelIso_zjets,"Z+jets"); |
| 195 |
> |
leg0->Draw(); |
| 196 |
> |
|
| 197 |
|
h_RelIso_all->Add(h_RelIso_ttbar,1.); |
| 198 |
|
h_RelIso_all->Add(h_RelIso_qcd,1.); |
| 199 |
|
h_RelIso_all->Add(h_RelIso_wjets,1.); |
| 200 |
|
h_RelIso_all->Add(h_RelIso_zjets,1.); |
| 201 |
|
h_RelIso_all->SetLineWidth(2); |
| 202 |
+ |
|
| 203 |
|
if( jbin == 4) |
| 204 |
|
h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 4 jets)"); |
| 205 |
|
else if ( jbin < 4 ) { |
| 208 |
|
title += " jets)"; |
| 209 |
|
h_RelIso_all->SetTitle(title); |
| 210 |
|
} |
| 211 |
< |
else {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;} |
| 211 |
> |
else { |
| 212 |
> |
cout<<"Don't mess around!!! jbin should be less than 5!"<<endl; |
| 213 |
> |
return; |
| 214 |
> |
} |
| 215 |
|
h_RelIso_all->SetMinimum(0); |
| 216 |
< |
// h_RelIso_all->Sumw2(); |
| 217 |
< |
// h_RelIso_all->Draw("e1"); |
| 218 |
< |
h_RelIso_all->Draw("same"); |
| 219 |
< |
|
| 216 |
> |
// Save for later use |
| 217 |
> |
TFile *file0 = new TFile("extpQCDbin4.root","RECREATE"); |
| 218 |
> |
file0->mkdir("histos"); |
| 219 |
> |
file0->cd("histos"); |
| 220 |
> |
h_RelIso_all->Write(); |
| 221 |
> |
h_RelIso_ttbar->Write(); |
| 222 |
> |
h_RelIso_qcd->Write(); |
| 223 |
> |
h_RelIso_wjets->Write(); |
| 224 |
> |
h_RelIso_zjets->Write(); |
| 225 |
> |
file0->Write(); |
| 226 |
> |
file0->Close(); |
| 227 |
|
// Fitting here: |
| 228 |
< |
FitBKG(fitfunc,h_RelIso_all,h_RelIso_qcd,cv1,AX0,AXf,BX0,BXf,Na,Nfac); |
| 192 |
< |
|
| 228 |
> |
FitBKG(fitfunc,h_RelIso_all,h_RelIso_qcd,AX0,AXf,BX0,BXf,Na); |
| 229 |
|
} |
| 230 |
|
|
| 231 |
< |
void qcd::FitBKG(TString s0,TH1D *h0,TH1D *h1,TCanvas *cv0, |
| 232 |
< |
double ax0,double axf,double bx0,double bxf,double na,double nfac) { |
| 233 |
< |
|
| 231 |
> |
void qcd::FitBKG(TString s0,TH1D *h0,TH1D *h1, |
| 232 |
> |
double ax0,double axf,double bx0,double bxf,double na) { |
| 233 |
> |
int niter = 1; |
| 234 |
> |
double pass[3] = {0.,0.,0.}; |
| 235 |
|
double ns, ns2; |
| 236 |
< |
ns = ns2 = 0; |
| 236 |
> |
ns = ns2 = 0.; |
| 237 |
> |
|
| 238 |
> |
TCanvas *cvf = new TCanvas("cvf","cvf",800,600); |
| 239 |
|
|
| 240 |
|
// Fit QCD in background region |
| 241 |
|
cout<<"Fitting with "<<s0<<" function!"<<endl; |
| 242 |
|
TH1D *h_all = (TH1D*)h0->Clone(); |
| 243 |
|
TH1D *h_qcd = (TH1D*)h1->Clone(); |
| 244 |
+ |
h_all->Draw(); |
| 245 |
+ |
|
| 246 |
+ |
cout << "\n>>>>> Iteration step: "<< niter << endl; |
| 247 |
+ |
h_all->Fit(s0,"0","",bx0,bxf); |
| 248 |
|
|
| 249 |
< |
h_all->Fit(s0,"","",bx0,bxf); |
| 249 |
> |
// Very first fit function |
| 250 |
> |
TF1 *f0 = (TF1*)h_all->GetFunction(s0); |
| 251 |
> |
Int_t npar = f0->GetNpar(); |
| 252 |
> |
Double_t chi2 = f0->GetChisquare(); |
| 253 |
> |
Int_t ndof = f0->GetNDF(); |
| 254 |
> |
pass[0] = chi2/ndof; |
| 255 |
> |
pass[1] = pass[0]; |
| 256 |
> |
cout << ">> Norm chi2 = " << pass[0] << endl; |
| 257 |
> |
|
| 258 |
> |
bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2); |
| 259 |
> |
bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2); |
| 260 |
|
|
| 208 |
– |
//Get a pointer to the fitted function |
| 209 |
– |
TF1 *f1 = (TF1*)h_all->GetFunction(s0); |
| 210 |
– |
Int_t npar = f1->GetNpar(); |
| 211 |
– |
// Double_t chi2 = f1->GetChisquare(); |
| 261 |
|
double * par0 = new double [npar]; |
| 262 |
|
double * parErr0 = new double [npar]; |
| 263 |
< |
printf("Function has %i parameters. Chisquare = %g\n", |
| 264 |
< |
npar, |
| 265 |
< |
f1->GetChisquare()); |
| 266 |
< |
for (Int_t i=0;i<npar;i++) { |
| 267 |
< |
printf("%s = %g +- %g\n", |
| 268 |
< |
f1->GetParName(i), |
| 269 |
< |
f1->GetParameter(i), |
| 270 |
< |
f1->GetParError(i) |
| 271 |
< |
); |
| 272 |
< |
par0[i] = f1->GetParameter(i); |
| 273 |
< |
parErr0[i] = f1->GetParError(i); |
| 274 |
< |
// cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl; |
| 263 |
> |
|
| 264 |
> |
double delta = pass[1]-pass[2]; |
| 265 |
> |
|
| 266 |
> |
while (niter <= 20 && abs(pass[1]-pass[2]) > 0.00001) { |
| 267 |
> |
if (niter > 1) |
| 268 |
> |
pass[1]=pass[2]; |
| 269 |
> |
if (delta >= 0){ |
| 270 |
> |
niter++; |
| 271 |
> |
cout << "\n>>>>> Iteration step: "<< niter << endl; |
| 272 |
> |
cout << ">> Pass[1] = " << pass[1] << endl; |
| 273 |
> |
|
| 274 |
> |
h_all->Fit(s0,"0","",bx0,bxf); |
| 275 |
> |
TF1 *fi = (TF1*)h_all->GetFunction(s0); |
| 276 |
> |
pass[2] = fi->GetChisquare()/fi->GetNDF(); |
| 277 |
> |
cout << ">> Norm chi2 = pass[2] = " << pass[2] << endl; |
| 278 |
> |
delta = pass[1]-pass[2]; |
| 279 |
> |
cout << ">> Delta = " << delta << endl; |
| 280 |
> |
|
| 281 |
> |
if (delta >= 0) { |
| 282 |
> |
printf("Function has %i parameters. Chisquare = %g\n", |
| 283 |
> |
npar, |
| 284 |
> |
fi->GetChisquare()); |
| 285 |
> |
for (Int_t i=0;i<npar;i++) { |
| 286 |
> |
printf("%s = %g +- %g\n", |
| 287 |
> |
fi->GetParName(i), |
| 288 |
> |
fi->GetParameter(i), |
| 289 |
> |
fi->GetParError(i) |
| 290 |
> |
); |
| 291 |
> |
par0[i] = fi->GetParameter(i); |
| 292 |
> |
parErr0[i] = fi->GetParError(i); |
| 293 |
> |
// cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl; |
| 294 |
> |
} |
| 295 |
> |
bx0 = fi->GetParameter(1) - wxmin*fi->GetParameter(2); |
| 296 |
> |
bxf = fi->GetParameter(1) + wxmax*fi->GetParameter(2); |
| 297 |
> |
cout << ">> bx0, bxf = " << bx0 << ", " << bxf << endl; |
| 298 |
> |
} |
| 299 |
> |
else { |
| 300 |
> |
cout << ">> Use previous fit!" << endl; |
| 301 |
> |
} |
| 302 |
> |
} |
| 303 |
> |
} |
| 304 |
> |
|
| 305 |
> |
cout << ">>> Final fitting resion (bx0, bxf) = " << bx0 << ", " << bxf << endl; |
| 306 |
> |
|
| 307 |
> |
// Get number of events within fitted region |
| 308 |
> |
TAxis *axis = h_all->GetXaxis(); |
| 309 |
> |
int bmin = axis->FindBin(bx0); |
| 310 |
> |
int bmax = axis->FindBin(bxf); |
| 311 |
> |
double nfac1 = h_all->Integral(bmin,bmax); |
| 312 |
> |
nfac1 -= (h_all->GetBinContent(bmin))*(bx0-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin); |
| 313 |
> |
nfac1 -= (h_all->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bxf)/axis->GetBinWidth(bmax); |
| 314 |
> |
cout << ">>> Final Nfac = " << nfac1 << endl; |
| 315 |
> |
|
| 316 |
> |
// //////////////////////// |
| 317 |
> |
// Histos and extrapolation |
| 318 |
> |
|
| 319 |
> |
TF1 *f1 = (TF1*)f0->Clone(); |
| 320 |
> |
f1->SetRange(bx0,bxf); |
| 321 |
> |
for (int i=0;i<npar;i++) { |
| 322 |
> |
cout<<" >> Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl; |
| 323 |
> |
f1->SetParameter(i,par0[i]); |
| 324 |
|
} |
| 325 |
+ |
|
| 326 |
|
f1->SetLineColor(2); |
| 327 |
|
h_qcd->SetLineColor(40); |
| 328 |
|
h_qcd->SetFillStyle(3018); |
| 329 |
|
h_qcd->SetFillColor(38); |
| 330 |
|
h_qcd->Draw("same"); |
| 331 |
|
f1->Draw("same"); |
| 332 |
< |
|
| 332 |
> |
|
| 333 |
> |
// Connect fitted and extrapolated region |
| 334 |
> |
TF1 *fin = (TF1*)f1->Clone(); |
| 335 |
> |
fin->SetRange(axf,bx0); |
| 336 |
> |
fin->SetLineStyle(2); |
| 337 |
> |
fin->SetLineColor(8); |
| 338 |
> |
fin->Draw("same"); |
| 339 |
> |
|
| 340 |
|
TF1 *f2 = (TF1*)f1->Clone(); |
| 341 |
|
f2->SetRange(ax0,axf); |
| 342 |
|
f2->SetLineColor(4); |
| 347 |
|
f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15); |
| 348 |
|
ns = f2->IntegralFast(np,x,w,ax0,axf); |
| 349 |
|
ns/=(f2->Integral(bx0,bxf)); |
| 350 |
< |
ns*=nfac; |
| 351 |
< |
cout<<"Ns by usual integral = "<<ns<<endl; |
| 350 |
> |
ns*=nfac1; |
| 351 |
> |
cout<<">>> Ns by usual integral = "<<ns<<endl; |
| 352 |
|
delete [] x; |
| 353 |
|
delete [] w; |
| 354 |
|
|
| 359 |
|
ig.SetFunction(wf); |
| 360 |
|
ns2 = ig.Integral(ax0,axf); |
| 361 |
|
ns2/=(g->Integral(bx0,bxf)); |
| 362 |
< |
ns2*=nfac; |
| 363 |
< |
cout<<"Ns by MathMore integral = "<<ns2<<endl; |
| 362 |
> |
ns2*=nfac1; |
| 363 |
> |
cout<<">>> Ns by MathMore integral = "<<ns2<<endl; |
| 364 |
|
|
| 365 |
< |
// Yet another estimation |
| 260 |
< |
double sum=0.,limitLow=ax0,limitHigh=axf,vx=0.,dx=0.; |
| 261 |
< |
dx=(limitHigh-limitLow)/nbin; |
| 262 |
< |
// cout<<"limitHigh = "<<limitHigh<<"; dx = "<<dx<<endl; |
| 263 |
< |
for(int i=1;i<nbin;i++){ |
| 264 |
< |
vx=h_all->GetXaxis()->GetBinCenter(i); |
| 265 |
< |
if(vx <= limitHigh){ |
| 266 |
< |
// cout<<"f2("<<vx<<")= "<<f2->Eval(vx)<<endl; |
| 267 |
< |
sum+=f2->Eval(vx); |
| 268 |
< |
} |
| 269 |
< |
} |
| 270 |
< |
cout<<"sum = "<<sum<<endl; |
| 271 |
< |
cv0->Update(); |
| 365 |
> |
cvf->Update(); |
| 366 |
|
|
| 367 |
|
TPaveStats *p1 = (TPaveStats*)h_all->GetListOfFunctions()->FindObject("stats"); |
| 368 |
|
p1->SetTextColor(kBlue); |
| 373 |
|
p1->Draw(); |
| 374 |
|
|
| 375 |
|
// Label histo |
| 376 |
< |
TLegend * leg1 = new TLegend(0.55,0.35,0.88,0.60); |
| 377 |
< |
leg1->SetHeader("Fit with "+fitfunc); |
| 376 |
> |
TLegend * leg1 = new TLegend(0.3,0.15,0.68,0.45); |
| 377 |
> |
// leg1->SetHeader("Fit with "+fitfunc); |
| 378 |
|
leg1->AddEntry(h_all,"All events (S+B)"); |
| 379 |
|
leg1->AddEntry(h_qcd,"QCD events"); |
| 380 |
< |
leg1->AddEntry(f1,"Fit to all events in control region"); |
| 381 |
< |
leg1->AddEntry(f2,"Extrapolation of QCD events in signal region"); |
| 380 |
> |
leg1->AddEntry(f1,"Fit of all events in control region"); |
| 381 |
> |
leg1->AddEntry(f2,"Extrapolation to signal region"); |
| 382 |
|
leg1->Draw(); |
| 383 |
|
// Print results |
| 384 |
< |
cout << " Observed Ns = " << ns; |
| 384 |
> |
cout << ">>>>> Observed Ns = " << ns; |
| 385 |
|
cout <<"\n"; |
| 386 |
< |
cout << " Expected Na = " << na << endl; |
| 386 |
> |
cout << ">>>>> Expected Na = " << na << endl; |
| 387 |
|
} |
