| 1 |
#if !defined(__CINT__) || defined(__MAKECINT__)
|
| 2 |
#include <TROOT.h> // access to gROOT, entry point to ROOT system
|
| 3 |
#include <TFile.h> // file handle class
|
| 4 |
#include <TTree.h> // class to access ntuples
|
| 5 |
#include <TChain.h>
|
| 6 |
#include <TCanvas.h> // class for drawing
|
| 7 |
#include <TH1F.h> // 1D histograms
|
| 8 |
#include <TH2F.h> // 2D histograms
|
| 9 |
#include <TBenchmark.h> // class to track macro running statistics
|
| 10 |
#include <TVector3.h> // 3D vector class
|
| 11 |
#include <TRegexp.h> // 3D vector class
|
| 12 |
#include <vector> // STL vector class
|
| 13 |
#include <iostream> // standard I/O
|
| 14 |
#include <iomanip> // functions to format standard I/O
|
| 15 |
#include <fstream> // functions for file I/O
|
| 16 |
#include <string> // C++ string class
|
| 17 |
#include <sstream> // class for parsing strings
|
| 18 |
#include "RooHistPdf.h"
|
| 19 |
#include "RooGlobalFunc.h"
|
| 20 |
#include "RooPlot.h"
|
| 21 |
#include "RooRealVar.h"
|
| 22 |
#include "RooDataHist.h"
|
| 23 |
#include "RooGenericPdf.h"
|
| 24 |
#include "RooDataSet.h"
|
| 25 |
#include "RooAddPdf.h"
|
| 26 |
#include "RooFormulaVar.h"
|
| 27 |
#include "RooCurve.h"
|
| 28 |
#include "RooBinning.h"
|
| 29 |
#include "Common/CPlot.hh" // helper class for plots
|
| 30 |
#include "Common/MitStyleRemix.hh" // style settings for drawing
|
| 31 |
#include "Common/MyTools.hh" // miscellaneous helper functions
|
| 32 |
#include "Common/CSample.hh" // helper class for organizing input ntuple files
|
| 33 |
|
| 34 |
#endif
|
| 35 |
using namespace RooFit;
|
| 36 |
|
| 37 |
TString scheme = "bdt-medium";
|
| 38 |
TString lumi = "2100";
|
| 39 |
TCanvas *c = new TCanvas("c","c");
|
| 40 |
vector<TString> fnamev,labelv;
|
| 41 |
|
| 42 |
void fitMet(TH1F *hmcZ, TH1F *hmcW, TH1F *hSSmet);
|
| 43 |
TH1F* get_FR_pred(TTree *tree, TString var, TString extra, int nbins, double xmin, double xmax ) {
|
| 44 |
TH1F* hFR = new TH1F("hFR","",nbins,xmin,xmax);
|
| 45 |
hFR->Sumw2();
|
| 46 |
tree->Draw(var+">>hFR","w*(w!=1)"+extra);
|
| 47 |
return hFR;
|
| 48 |
}
|
| 49 |
|
| 50 |
TH1F* get_SS(TTree *tree, TString var, TString extra, int nbins, double xmin, double xmax) {
|
| 51 |
TH1F* hSS = new TH1F("hSS","",nbins,xmin,xmax);
|
| 52 |
hSS->Sumw2();
|
| 53 |
tree->Draw(var+">>hSS","(w==1)"+extra);
|
| 54 |
return hSS;
|
| 55 |
}
|
| 56 |
|
| 57 |
TH1F* get_MC(TTree *tree, TString var, TString extra, int nbins, double xmin, double xmax) {
|
| 58 |
TH1F *hMC = new TH1F("hMC","",nbins,xmin,xmax);
|
| 59 |
hMC->Sumw2();
|
| 60 |
tree->Draw(var+">>hMC",lumi+"*npuw*xswgt*(w==1)"+extra);
|
| 61 |
return hMC;
|
| 62 |
}
|
| 63 |
|
| 64 |
void plot_fake_closure()
|
| 65 |
{
|
| 66 |
CPlot::sOutDir = "/home/dkralph/public_html/4lplots";
|
| 67 |
|
| 68 |
map<TString,TTree*> treemap;
|
| 69 |
map<TString,TString> labelmap;
|
| 70 |
|
| 71 |
TChain *wdata_chain = new TChain("wdata_chain","wdata_chain");
|
| 72 |
wdata_chain->AddFile("root/ntuples/W-bdt-medium-test.root",-1,"passtuple");
|
| 73 |
treemap["wdata"] = wdata_chain;
|
| 74 |
|
| 75 |
TChain *zmc_chain = new TChain("zmc_chain","zmc_chain");
|
| 76 |
zmc_chain->AddFile("root/ntuples/W-bdt-medium-zjets-test.root",-1,"passtuple");
|
| 77 |
treemap["zmc"] = zmc_chain;
|
| 78 |
|
| 79 |
TChain *wmc_chain = new TChain("wmc_chain","wmc_chain");
|
| 80 |
wmc_chain->AddFile("root/ntuples/W-bdt-medium-wjets-test.root",-1,"passtuple");
|
| 81 |
treemap["wmc"] = wmc_chain;
|
| 82 |
|
| 83 |
double xmin,xmax;
|
| 84 |
int nbins;
|
| 85 |
TString var,cut,extra,ylabel;
|
| 86 |
|
| 87 |
//----------------------------------------------------------------------------------------
|
| 88 |
// EMu W selection
|
| 89 |
cut = "*(pttight>25 && pfisotight/pttight<.025)";
|
| 90 |
ylabel = "events";
|
| 91 |
|
| 92 |
var = "mZ1";
|
| 93 |
xmin = 0; xmax = 300; nbins = 20;
|
| 94 |
extra = cut + "*(met>30)";
|
| 95 |
CPlot plotmZ1W(var+"W","","#bf{m_{ll} [GeV]}",ylabel);
|
| 96 |
TH1F *hSSmZ1 = get_SS(treemap["wdata"],var,extra,nbins,xmin,xmax);
|
| 97 |
plotmZ1W.AddHist1D(hSSmZ1,"Obs.","E");
|
| 98 |
TH1F *hFRmZ1 = get_FR_pred(treemap["wdata"],var,extra,nbins,xmin,xmax);
|
| 99 |
plotmZ1W.AddToStack(hFRmZ1,"FR predic",kCyan-6);
|
| 100 |
plotmZ1W.Draw(c,kTRUE,"png");
|
| 101 |
|
| 102 |
extra = cut + "*(met>30 && pt2<20)";
|
| 103 |
CPlot plotmZ1W_lopt(var+"W_lopt","","#bf{m_{ll} [GeV]}",ylabel);
|
| 104 |
TH1F *hSSmZ1_lopt = get_SS(treemap["wdata"],var,extra,nbins,xmin,xmax);
|
| 105 |
plotmZ1W_lopt.AddHist1D(hSSmZ1_lopt,"Obs.","E");
|
| 106 |
TH1F *hFRmZ1_lopt = get_FR_pred(treemap["wdata"],var,extra,nbins,xmin,xmax);
|
| 107 |
plotmZ1W_lopt.AddToStack(hFRmZ1_lopt,"FR predic",kCyan-6);
|
| 108 |
plotmZ1W_lopt.Draw(c,kTRUE,"png");
|
| 109 |
|
| 110 |
extra = cut + "*(met>30 && pt2>20)";
|
| 111 |
CPlot plotmZ1W_hipt(var+"W_hipt","","#bf{m_{ll} [GeV]}",ylabel);
|
| 112 |
TH1F *hSSmZ1_hipt = get_SS(treemap["wdata"],var,extra,nbins,xmin,xmax);
|
| 113 |
plotmZ1W_hipt.AddHist1D(hSSmZ1_hipt,"Obs.","E");
|
| 114 |
TH1F *hFRmZ1_hipt = get_FR_pred(treemap["wdata"],var,extra,nbins,xmin,xmax);
|
| 115 |
plotmZ1W_hipt.AddToStack(hFRmZ1_hipt,"FR predic",kCyan-6);
|
| 116 |
plotmZ1W_hipt.Draw(c,kTRUE,"png");
|
| 117 |
|
| 118 |
var = "pt2";
|
| 119 |
xmin = 0; xmax = 80; nbins = 30;
|
| 120 |
extra = cut;
|
| 121 |
CPlot plotpt2W(var+"W","","#bf{loose p_{T} [GeV]}",ylabel);
|
| 122 |
TH1F *pt2_SS = get_SS(treemap["wdata"],var,extra,nbins,xmin,xmax);
|
| 123 |
plotpt2W.AddHist1D(pt2_SS,"Obs.","E");
|
| 124 |
TH1F *pt2_FR = get_FR_pred(treemap["wdata"],var,extra,nbins,xmin,xmax);
|
| 125 |
plotpt2W.AddToStack(pt2_FR,"FR predic",kCyan-6);
|
| 126 |
plotpt2W.Draw(c,kTRUE,"png");
|
| 127 |
for(int i=0; i<pt2_SS->GetNbinsX()+2; i++) {
|
| 128 |
cout << pt2_SS->GetBinContent(i) << "\t" << pt2_FR->GetBinContent(i) << "\t" << (pt2_SS->GetBinContent(i) - pt2_FR->GetBinContent(i))/(0.5*(pt2_SS->GetBinContent(i) + pt2_FR->GetBinContent(i))) << endl;
|
| 129 |
}
|
| 130 |
|
| 131 |
var = "abs(eta2)";
|
| 132 |
xmin = 0; xmax = 2.6; nbins = 30;
|
| 133 |
extra = cut;
|
| 134 |
CPlot ploteta2W("eta2W","","#bf{loose p_{T} [GeV]}",ylabel);
|
| 135 |
TH1F *htmp = get_SS(treemap["wdata"],var,extra,nbins,xmin,xmax);
|
| 136 |
ploteta2W.AddHist1D(htmp,"Obs.","E");
|
| 137 |
ploteta2W.AddToStack(get_FR_pred(treemap["wdata"],var,extra,nbins,xmin,xmax),"FR predic",kCyan-6);
|
| 138 |
ploteta2W.SetYRange(0,1.4*htmp->GetMaximum());
|
| 139 |
ploteta2W.Draw(c,kTRUE,"png");
|
| 140 |
|
| 141 |
var = "mt";
|
| 142 |
xmin = 0; xmax = 120; nbins = 30;
|
| 143 |
extra = cut + "*(met>30)";
|
| 144 |
CPlot plotMtW(var+"W","","#bf{m_{T} [GeV]}",ylabel);
|
| 145 |
plotMtW.AddHist1D(get_SS(treemap["wdata"],var,extra,nbins,xmin,xmax),"Obs.","E");
|
| 146 |
plotMtW.AddToStack(get_FR_pred(treemap["wdata"],var,extra,nbins,xmin,xmax),"FR predic",kCyan-6);
|
| 147 |
plotMtW.Draw(c,kTRUE,"png");
|
| 148 |
|
| 149 |
//----------------------------------------------------------------------------------------
|
| 150 |
// met fit
|
| 151 |
var = "met";
|
| 152 |
xmin = 0; xmax = 100; nbins = 25;
|
| 153 |
extra = cut;
|
| 154 |
TH1F *hSSmet = get_SS(treemap["wdata"],var,extra,nbins,xmin,xmax);
|
| 155 |
TH1F *hmcZ = get_MC(treemap["zmc"],var,extra,nbins,xmin,xmax); hmcZ->Scale(1./6.5);
|
| 156 |
TH1F *hmcW = get_MC(treemap["wmc"],var,extra,nbins,xmin,xmax); hmcW->Scale(1./6.5);
|
| 157 |
|
| 158 |
fitMet(hmcZ,hmcW,hSSmet);
|
| 159 |
|
| 160 |
double integral=0,err=0;
|
| 161 |
vector<double> integrals,errs;
|
| 162 |
// inclusive
|
| 163 |
integral = hSSmZ1->IntegralAndError(1,hSSmZ1->GetNbinsX(),err); integrals.push_back(integral); errs.push_back(err);
|
| 164 |
integral = hFRmZ1->IntegralAndError(1,hFRmZ1->GetNbinsX(),err); integrals.push_back(integral); errs.push_back(err);
|
| 165 |
// for(int i=0; i<hSSmZ1->GetNbinsX()+2; i++) {
|
| 166 |
// cout << hSSmZ1->GetBinContent(i) << "\t" << hFRmZ1->GetBinContent(i) << "\t" << (hSSmZ1->GetBinContent(i) - hFRmZ1->GetBinContent(i))/(0.5*(hSSmZ1->GetBinContent(i) + hFRmZ1->GetBinContent(i))) << endl;
|
| 167 |
// }
|
| 168 |
// low pt
|
| 169 |
integral = hSSmZ1_lopt->IntegralAndError(1,hSSmZ1_lopt->GetNbinsX(),err); integrals.push_back(integral); errs.push_back(err);
|
| 170 |
integral = hFRmZ1_lopt->IntegralAndError(1,hFRmZ1_lopt->GetNbinsX(),err); integrals.push_back(integral); errs.push_back(err);
|
| 171 |
// high pt
|
| 172 |
integral = hSSmZ1_hipt->IntegralAndError(1,hSSmZ1_hipt->GetNbinsX(),err); integrals.push_back(integral); errs.push_back(err);
|
| 173 |
integral = hFRmZ1_hipt->IntegralAndError(1,hFRmZ1_hipt->GetNbinsX(),err); integrals.push_back(integral); errs.push_back(err);
|
| 174 |
|
| 175 |
cout << "Inclusive: " << endl;
|
| 176 |
cout << "\tSS: " << fixed << setprecision(1) << integrals[0] << "\t+/- " << fixed << setprecision(1) << setw(5) << errs[0] << "\t (events: " << hSSmZ1->GetEntries() << ")" << endl;
|
| 177 |
cout << "\tFR: " << fixed << setprecision(1) << integrals[1] << "\t+/- " << fixed << setprecision(1) << setw(5) << errs[1] << "\t (events: " << hFRmZ1->GetEntries() << ")\t";
|
| 178 |
cout << setprecision(3) << 2*(integrals[1]-integrals[0])/(integrals[0]+integrals[1]) << endl;
|
| 179 |
cout << "Pt < 20: " << endl;
|
| 180 |
cout << "\tSS: " << fixed << setprecision(1) << integrals[2] << "\t+/- " << fixed << setprecision(1) << setw(5) << errs[2] << "\t (events: " << hSSmZ1_lopt->GetEntries() << ")" << endl;
|
| 181 |
cout << "\tFR: " << fixed << setprecision(1) << integrals[3] << "\t+/- " << fixed << setprecision(1) << setw(5) << errs[3] << "\t (events: " << hFRmZ1_lopt->GetEntries() << ")\t";
|
| 182 |
cout << setprecision(3) << 2*(integrals[3]-integrals[2])/(integrals[3]+integrals[2]) << endl;
|
| 183 |
cout << "Pt > 20: " << endl;
|
| 184 |
cout << "\tSS: " << fixed << setprecision(1) << integrals[4] << "\t+/- " << fixed << setprecision(1) << setw(5) << errs[4] << "\t (events: " << hSSmZ1_hipt->GetEntries() << ")" << endl;
|
| 185 |
cout << "\tFR: " << fixed << setprecision(1) << integrals[5] << "\t+/- " << fixed << setprecision(1) << setw(5) << errs[5] << "\t (events: " << hFRmZ1_hipt->GetEntries() << ")\t";
|
| 186 |
cout << setprecision(3) << 2*(integrals[5]-integrals[4])/(integrals[5]+integrals[4]) << endl;
|
| 187 |
}
|
| 188 |
//----------------------------------------------------------------------------------------
|
| 189 |
void fitMet(TH1F *hmcZ, TH1F *hmcW, TH1F *hSSmet)
|
| 190 |
{
|
| 191 |
RooRealVar met("met","MET [GeV]",0,100);
|
| 192 |
RooBinning metbins(25,0,100);
|
| 193 |
met.setBinning(metbins);
|
| 194 |
// Z mc pdf
|
| 195 |
RooDataHist hZ("hZ","hZ",RooArgList(met),hmcZ);
|
| 196 |
RooHistPdf pZ("pZ","pZ",RooArgSet(met),hZ);
|
| 197 |
// W mc pdf
|
| 198 |
hmcW->Smooth(2);
|
| 199 |
RooDataHist hW("hW","hW",RooArgList(met),hmcW);
|
| 200 |
RooHistPdf pW("pW","pW",RooArgSet(met),hW,3);
|
| 201 |
// rayleigh pdf
|
| 202 |
RooRealVar sig("sig","sig",3,1,15);
|
| 203 |
RooRealVar a1("a1","a1",0.1,0,1);
|
| 204 |
RooGenericPdf rayl("rayl","rayl","met*exp(-(met)*(met)/(sig*sig +2*sig*a1*(met) + a1*a1*(met)*(met)))",RooArgList(met,sig,a1));
|
| 205 |
|
| 206 |
RooRealVar nWandZ("nWandZ","nWandZ",10000,0,11000);//81+611,1000);
|
| 207 |
RooRealVar fZ("fZ","fZ",84./628);
|
| 208 |
RooRealVar nRayl("nRayl","nRayl",300,0,10000);
|
| 209 |
RooAddPdf mMC("mMC","mMC",RooArgList(pZ,pW),RooArgList(fZ));
|
| 210 |
RooAddPdf model("model","model",RooArgList(mMC,rayl),RooArgList(nWandZ,nRayl));
|
| 211 |
|
| 212 |
RooDataHist obsMet("obsMet","obsMet",RooArgList(met),hSSmet);
|
| 213 |
model.fitTo(obsMet);
|
| 214 |
|
| 215 |
CPlot plotMetFit("metfitW","","MET [GeV]","Events");
|
| 216 |
TH1F *rayl_hist = (TH1F*)rayl.createHistogram("met",met,Binning(metbins));
|
| 217 |
rayl_hist->Scale(nRayl.getVal()/rayl_hist->Integral());
|
| 218 |
plotMetFit.AddToStack(rayl_hist,"jet (Rayleigh)",kBlue);
|
| 219 |
TH1F *z_hist = (TH1F*)pZ.createHistogram("met",met,Binning(metbins));
|
| 220 |
z_hist->Scale(fZ.getVal()*nWandZ.getVal()/z_hist->Integral());
|
| 221 |
plotMetFit.AddToStack(z_hist,"Z+jets",kGray);
|
| 222 |
TH1F *w_hist = (TH1F*)pW.createHistogram("met",met,Binning(metbins));
|
| 223 |
w_hist->Scale( (nWandZ.getVal()*(1-fZ.getVal()))/w_hist->Integral());
|
| 224 |
plotMetFit.AddToStack(w_hist,"W+jets",kCyan-6);
|
| 225 |
TH1F *h_obs = (TH1F*)obsMet.createHistogram("met",met,Binning(metbins));
|
| 226 |
plotMetFit.AddHist1D(h_obs,"Obs.");
|
| 227 |
plotMetFit.AddLine(30,0,30,850,kRed,kDashed);
|
| 228 |
plotMetFit.Draw(c,kTRUE,"png");
|
| 229 |
|
| 230 |
|
| 231 |
cout << "\nnW in mc: " << hmcW->Integral() << endl;
|
| 232 |
cout << "nZ in mc: " << hmcZ->Integral() << endl << endl;
|
| 233 |
}
|
| 234 |
|
