ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Jeng/tools/ikstest.cc
Revision: 1.5
Committed: Wed Sep 29 23:40:33 2010 UTC (14 years, 7 months ago) by jengbou
Content type: text/plain
Branch: MAIN
Changes since 1.4: +115 -80 lines
Log Message:
update

File Contents

# User Rev Content
1 jengbou 1.1 #include "ikstest.h"
2 jengbou 1.5 #include "TopStyle/tdrstyle.C"
3 jengbou 1.1 #include "TopStyle/CMSTopStyle.cc"
4    
5     //=================================
6     // Program to run IKS test
7     // * Input directories:
8     // * Data: skimmed_Data_x.xxpb-1
9     // => default selection : Data_<suffix>.root
10     // => n-1 selection : Data_<suffix>_<invName>.root
11     // * MC : skimmed_MC {QCD,WJets,TTbar...}
12     // => default selection : QCD_<suffix>.root
13     // => n-1 selection : Data_<suffix>_<invName>.root
14     // * Upmost output dir specified by <desDir>
15     // subdirectory created according to <useInv> and <realData>
16     // => <desDir'><xyz><suffix>.pdf (.txt)
17     // * e.g. to do KS test on data with MC QCD shape in signal region:
18     // useInv=false; realData=true
19     //=================================
20    
21     using namespace std;
22    
23     //define the constants: 2.88/pb
24     const Double_t weight_[5] = {0.0524313, //QCD
25     0.0089567, //WJets
26     0.000306, //TTbar
27     0.0080911, //ZJets
28     0.000114 //STtch
29     };
30     const Double_t procQCD = 1.46;
31     const Double_t procWJets = 1.03;
32     const Double_t procTTbar = 1.;
33     const Double_t procZJets = 1.;
34     const Double_t procSTtch = 1.;
35    
36     // Output directory
37 jengbou 1.3 TString baseDir = "Results_2.88pb-1_";
38 jengbou 1.5 //TString jbinNames[] = {"Incl0","Incl1","Incl2","Incl3","Incl4","Excl0","Excl1","Excl2","Excl3"};
39     TString jbinNames[] = {"Incl0"};
40     TString postDir = "_TEST0929/";
41 jengbou 1.1 // User defined parameters
42 jengbou 1.3 bool realData = true;
43     bool useInv = true; // whether to use n-1 QCD template
44     bool plotInvSelShape = false; // whether to use shape from inverse selection for dataKS plots. Available when useInv is true
45 jengbou 1.5 bool useIKSopt = true;
46     double yLowBound = 1.e-36;
47 jengbou 1.1 // Ntuples to use
48 jengbou 1.5 TString selNames[] = {"Sel5"}; // Suffix of selection
49     // TString invNames[] = {"METlt10RelIsogt0p1"};
50     // TString invLabels[] = {"#slash{E}_{T} < 10, RelIso > 0.1"};
51     TString invNames[] = {"RelIsogt0p1","D0gt0p02"};
52     TString invLabels[] = {"RelIso > 0.1","d_{0} > 0.02"};
53 jengbou 1.3 bool debug_ = false;
54 jengbou 1.4 bool UseOverflow = true;
55     bool exclusive_ = false;
56 jengbou 1.5 bool overwriteAll = false;
57    
58 jengbou 1.1 //=================================
59     // Main program
60     //=================================
61     void ikstest() {
62     CMSTopStyle style;
63 jengbou 1.5 setTDRStyle();
64     style.setupICHEPv1();
65 jengbou 1.1 TLatex *latex = new TLatex();
66     latex->SetNDC();
67    
68 jengbou 1.4 int n_sels = sizeof(selNames)/sizeof(selNames[0]);
69     int njbins_ = sizeof(jbinNames)/sizeof(jbinNames[0]);
70     int size_ninv = (useInv ? sizeof(invNames)/sizeof(invNames[0]) : 1);
71    
72     // if (debug_) {
73     cout << "size of selNames = " << n_sels << endl;
74     cout << "size of jbinNames = " << njbins_ << endl;
75     cout << "size of invNames = " << size_ninv << endl;
76     // }
77    
78     for (int i = 0; i < njbins_; ++i){
79 jengbou 1.3 struct stat stbDir;
80 jengbou 1.4 TString tmpDir = baseDir + jbinNames[i] + postDir;
81 jengbou 1.3 if (stat(tmpDir,&stbDir)){
82     cout << "Creating folder : " << tmpDir << endl;
83     if (mkdir(tmpDir,0755) == -1){
84     std::cerr << "Error creating folder" << endl;
85     return;
86     }
87 jengbou 1.1 }
88     }
89    
90     for (int ninv = 0; ninv < size_ninv; ++ninv) {
91     TString invName = invNames[ninv];
92     for (int nsel = 0; nsel < n_sels; ++nsel) {
93 jengbou 1.3 ostringstream ssel_;
94     ssel_ << nsel;
95     TString suffix0_ = selNames[nsel];
96 jengbou 1.4
97 jengbou 1.1 //open the files with histograms
98     map<string,TFile*> mfile;
99 jengbou 1.3 map<TString,TCanvas*> cvs; // map of usual histogram
100 jengbou 1.1 mfile.clear();
101 jengbou 1.3 cvs.clear();
102     mfile["Data"] = TFile::Open(TString("skimmed_Data_2.88pb-1_M3/Data_"+suffix0_+".root"));
103 jengbou 1.1 // n-1 cuts
104     if (useInv) {
105     if (realData)
106 jengbou 1.3 mfile["InvSel"] = TFile::Open(TString("skimmed_Data_2.88pb-1_M3/Data_"+suffix0_+"_"+invName+".root"));
107 jengbou 1.1 else
108 jengbou 1.3 mfile["InvSel"] = TFile::Open(TString("skimmed_MC_M3/v5/QCD_"+suffix0_+"_"+invName+".root"));
109 jengbou 1.1 }
110     // RefSel MC
111 jengbou 1.3 mfile["0"] = TFile::Open(TString("skimmed_MC_M3/v5/QCD_"+suffix0_+".root"));
112     mfile["1"] = TFile::Open(TString("skimmed_MC_M3/v5/WJets_"+suffix0_+".root"));
113     mfile["2"] = TFile::Open(TString("skimmed_MC_M3/v5/TTbar_"+suffix0_+".root"));
114     mfile["3"] = TFile::Open(TString("skimmed_MC_M3/v5/ZJets_"+suffix0_+".root"));
115     mfile["4"] = TFile::Open(TString("skimmed_MC_M3/v5/STtch_"+suffix0_+".root"));
116    
117 jengbou 1.4 for (int njbin = 0; njbin < njbins_; ++njbin) {
118     exclusive_ = false;
119     if (jbinNames[njbin].Contains("Excl")) exclusive_ = true;
120     char nthjbin_ = jbinNames[njbin](jbinNames[njbin].Length()-1);
121    
122 jengbou 1.3 TString desDir;
123 jengbou 1.4 TString tmpDir = baseDir + jbinNames[njbin] + postDir;
124 jengbou 1.3 // create output directory
125     if (!useInv) {
126     if (!realData) desDir = tmpDir + "MC/";
127     else desDir = tmpDir + "Data_MC/";
128     } else {
129     if (!realData) desDir = tmpDir + "MC_"+invName+"/";
130     else desDir = tmpDir + "Data_"+invName+"/";
131     }
132    
133     if (nsel == 0) {
134     struct stat stDir;
135     if (!stat(desDir,&stDir)){
136     cout << "\nOutput folder \"" << desDir
137     << "\" exists! Continue to overwrite it? (enter to continue; 'q' for quit)" << endl;
138     char incmd;
139 jengbou 1.4 if (overwriteAll) incmd='a';
140     else cin.get(incmd);
141 jengbou 1.3 if (incmd == 'q') return;
142     } else {
143     cout << "Creating folder : " << desDir << endl;
144     if (mkdir(desDir,0755) == -1){
145     std::cerr << "Error creating folder" << endl;
146     return;
147     }
148     }
149 jengbou 1.1 }
150    
151 jengbou 1.4 TString suffix1_ = jbinNames[njbin] + "_" + selNames[nsel];
152 jengbou 1.1
153 jengbou 1.3 ofstream outprint(TString(desDir+"Results_"+suffix1_+".txt"));
154 jengbou 1.1
155 jengbou 1.3 //define histograms and related parameters
156 jengbou 1.5 string histoName[] = {"h_mu_pt","h_met","h_mt"};
157     // string histoName[] = {"h_mu_pt","h_mu_eta","h_Ht"};
158     const int nvars_ = sizeof(histoName)/sizeof(histoName[0]);
159     if (nvars_ == 1) useIKSopt = false;
160     cout << "===> Fit " << nvars_ << " variable(s): " << endl;
161    
162     string histoTitle[nvars_];
163     string histoLabelX[nvars_];
164 jengbou 1.4
165     //Histo binning
166 jengbou 1.5 Int_t xbins[nvars_];
167     Double_t xlow[nvars_];
168     Double_t xhigh[nvars_];
169 jengbou 1.3 string sample[5] = {"QCD","WJets","TTbar","ZJets","STtch"};
170    
171 jengbou 1.5 TH1F* hMC_[5][nvars_]; // MC histograms
172     TH1F* hData_[nvars_]; // Data or mix MC
173     TH1F* hQCD_NEW[nvars_]; // InvSel QCD shape
174     TH1F* hKSres_[nvars_];
175     TH1F* hKSvalues_[nvars_+1];
176     TH1F* hQCD_KS[nvars_];
177     TH1F* hWJets_KS[nvars_];
178 jengbou 1.3
179 jengbou 1.5 if (useIKSopt) hKSvalues_[nvars_] = new TH1F(TString("productKSvalues_"+suffix1_),"",2./stepsize, stepsize, 2.+stepsize);
180 jengbou 1.4
181     string suffixj_ = "_";
182     if (exclusive_) suffixj_ += "excl";
183 jengbou 1.3 //load the histograms from the root files
184 jengbou 1.5 for (int vi = 0; vi < nvars_; ++vi) {// nvars_ variables
185 jengbou 1.3 //cout << "file[" << vi << "] : " << endl;
186     ostringstream svi_;
187     svi_ << vi;
188    
189     if (!useInv) {//use QCD MC sample
190 jengbou 1.4 hQCD_NEW[vi] = (TH1F*) mfile["0"]->Get(TString(histoName[vi]+suffixj_+nthjbin_))->Clone();
191 jengbou 1.3 hQCD_NEW[vi] -> Scale(weight_[0]);
192 jengbou 1.4 hQCD_NEW[vi] -> SetName(TString("InvSel_"+histoName[vi]+suffixj_+svi_.str()+"_"+nthjbin_));
193 jengbou 1.3 }
194     else {
195 jengbou 1.4 hQCD_NEW[vi] = (TH1F*) mfile["InvSel"]->Get(TString(histoName[vi]+suffixj_+nthjbin_))->Clone();
196 jengbou 1.3 if (!realData) hQCD_NEW[vi] -> Scale(weight_[0]);
197 jengbou 1.4 hQCD_NEW[vi] -> SetName(TString("InvSel_"+histoName[vi]+suffixj_+svi_.str()+"_"+nthjbin_));
198 jengbou 1.3 }
199 jengbou 1.5
200     xbins[vi] = hQCD_NEW[vi]->GetNbinsX();
201     xlow[vi] = hQCD_NEW[vi]->GetXaxis()->GetXmin();
202     xhigh[vi] = hQCD_NEW[vi]->GetXaxis()->GetXmax();
203     histoTitle[vi] = hQCD_NEW[vi]->GetTitle();
204     histoLabelX[vi] = hQCD_NEW[vi]->GetXaxis()->GetTitle();
205     cout << histoLabelX[vi] << ": nbins = " << xbins[vi] << "; x = (" << xlow[vi] << ", " << xhigh[vi] << ")" << endl;
206    
207 jengbou 1.3 if (debug_) cout << "hQCD_NEW[" << vi << "] @ " << hQCD_NEW[vi] << endl;
208 jengbou 1.5 if (UseOverflow) {
209 jengbou 1.4 hQCD_NEW[vi]->SetBinContent(xbins[vi],
210     hQCD_NEW[vi]->GetBinContent(xbins[vi]+1)+
211     hQCD_NEW[vi]->GetBinContent(xbins[vi]));
212 jengbou 1.5 double err_ = sqrt(hQCD_NEW[vi]->GetBinError(xbins[vi]+1)*hQCD_NEW[vi]->GetBinError(xbins[vi]+1)+
213     hQCD_NEW[vi]->GetBinError(xbins[vi])*hQCD_NEW[vi]->GetBinError(xbins[vi]));
214     hQCD_NEW[vi]->SetBinError(xbins[vi],err_);
215     }
216     string nameNewHisto = "mix_"+histoName[vi]+suffixj_+nthjbin_;
217     string nameNewHistoSFKS = "finalSF_"+histoName[vi]+suffixj_+nthjbin_;
218     string nameNewHistoKSvalues = "KSvalues_"+histoLabelX[vi]+suffixj_+nthjbin_;
219     hData_[vi] = new TH1F(nameNewHisto.c_str(),"",xbins[vi],xlow[vi],xhigh[vi]);
220     hKSres_[vi] = new TH1F(nameNewHistoSFKS.c_str(),"",xbins[vi],xlow[vi],xhigh[vi]);
221     hKSvalues_[vi] = new TH1F(nameNewHistoKSvalues.c_str(),"",2./stepsize, stepsize, 2.+stepsize);
222 jengbou 1.4 hData_[vi]->Sumw2();
223     hKSres_[vi]->Sumw2();
224     hKSvalues_[vi]->Sumw2();
225 jengbou 1.3 }
226 jengbou 1.5 //return;
227 jengbou 1.3
228     for (int n = 0; n < 5; ++n) {// MC samples
229 jengbou 1.5 for (int ihisto = 0; ihisto < nvars_; ihisto++) {// nvars_ variables
230 jengbou 1.3 //cout << "Variable[" << ihisto << "]" << endl;
231 jengbou 1.4 string histo_name = histoName[ihisto]+suffixj_+nthjbin_+"_"+sample[n];
232 jengbou 1.3 ostringstream ss_;
233     ss_ << n;
234 jengbou 1.4 //cout << TString(histoName[ihisto]+suffixj_+nthjbin_) << endl;
235     hMC_[n][ihisto] = (TH1F*) mfile[ss_.str()]->Get(TString(histoName[ihisto]+suffixj_+nthjbin_))->Clone();
236 jengbou 1.5 if (UseOverflow) {
237 jengbou 1.4 hMC_[n][ihisto]->SetBinContent(xbins[ihisto],
238     hMC_[n][ihisto]->GetBinContent(xbins[ihisto]+1)+
239     hMC_[n][ihisto]->GetBinContent(xbins[ihisto]));
240 jengbou 1.5 double err_ = sqrt(hMC_[n][ihisto]->GetBinError(xbins[ihisto]+1)*hMC_[n][ihisto]->GetBinError(xbins[ihisto]+1)+
241     hMC_[n][ihisto]->GetBinError(xbins[ihisto])*hMC_[n][ihisto]->GetBinError(xbins[ihisto]));
242     hMC_[n][ihisto]->SetBinError(xbins[ihisto],err_);
243     }
244 jengbou 1.3 if (debug_) {
245     cout << "File[" << n << "] @" << mfile[ss_.str()]
246 jengbou 1.4 << "; histo[" << ihisto << "] @ " << mfile[ss_.str()]->Get(TString(histoName[ihisto]+suffixj_+nthjbin_))
247 jengbou 1.3 <<"; hMC_[" << n << "][" << ihisto << "] raw evts = "
248     << setw(12) << hMC_[n][ihisto]->Integral();
249     }
250 jengbou 1.4
251 jengbou 1.3 hMC_[n][ihisto] -> Scale(weight_[n]);
252 jengbou 1.4 hMC_[n][ihisto] -> SetName(TString("MC_"+sample[n]+"_"+histoName[ihisto]+suffixj_+nthjbin_));
253 jengbou 1.3 if (debug_) cout << "; weighted num evts = " << setw(8) << hMC_[n][ihisto]->Integral() << endl;
254     }
255 jengbou 1.1 }
256    
257 jengbou 1.3 //create the mixed samples = "data"
258     TString cvsName0 = "Data";
259     if (useInv) {
260     cvsName0 += "_";
261     cvsName0 += invName;
262     }
263 jengbou 1.4 cvsName0 += TString("_" + ssel_.str() + nthjbin_);
264 jengbou 1.5 cvs[cvsName0] = new TCanvas(cvsName0,"Data distributions",1400,700);
265     cvs[cvsName0]->Divide(nvars_,1);
266     for (int i = 0; i < nvars_; i++) {
267 jengbou 1.3 cvs[cvsName0]->cd(i+1);
268     if (!realData) {
269     hData_[i] -> Add(hMC_[0][i],hMC_[1][i], procQCD,procWJets);
270     hData_[i] -> Add(hMC_[2][i], procTTbar);
271     hData_[i] -> Add(hMC_[3][i], procZJets);
272     hData_[i] -> Add(hMC_[4][i], procSTtch);
273     }
274     else {
275 jengbou 1.4 TH1F *htmp = (TH1F*) mfile["Data"]->Get(TString(histoName[i]+suffixj_+nthjbin_));
276 jengbou 1.3 hData_[i] -> Add(htmp,1.);
277 jengbou 1.5 if (UseOverflow) {
278 jengbou 1.4 hData_[i]->SetBinContent(xbins[i],
279     hData_[i]->GetBinContent(xbins[i]+1)+
280     hData_[i]->GetBinContent(xbins[i]));
281 jengbou 1.5 double err_ = sqrt(hData_[i]->GetBinError(xbins[i]+1)*hData_[i]->GetBinError(xbins[i]+1)+
282     hData_[i]->GetBinError(xbins[i])*hData_[i]->GetBinError(xbins[i]));
283     hData_[i]->SetBinError(xbins[i],err_);
284     }
285 jengbou 1.3 }
286     hData_[i]->GetXaxis()->SetTitle(histoLabelX[i].c_str());
287     hData_[i]->GetYaxis()->SetTitle("Events");
288     if (i == 0) // mu pT plot
289     hData_[0]->GetXaxis()->SetRangeUser(20.,xhigh[0]);
290     hData_[i]->DrawClone("hist");
291     }
292 jengbou 1.4 cvs[cvsName0]->SaveAs(TString(desDir+"Data_distributions_"+selNames[nsel]+".pdf"));
293     cvs[cvsName0]->SaveAs(TString(desDir+"Data_distributions_"+selNames[nsel]+".png"));
294 jengbou 1.3
295     //Calculate num of events for each sample
296     vector<double> vNev_;
297 jengbou 1.4 double errNevt[8];
298 jengbou 1.5 double NevData = hData_[nvars_-1]->IntegralAndError(1,xbins[nvars_-1],errNevt[0]);
299     double corr_NevQCD = hMC_[0][nvars_-1]->IntegralAndError(1,xbins[nvars_-1],errNevt[1]);
300     double corr_NevQCD_NEW = hQCD_NEW[nvars_-1]->IntegralAndError(1,xbins[nvars_-1],errNevt[2]);
301     double corr_NevWJets = hMC_[1][nvars_-1]->IntegralAndError(1,xbins[nvars_-1],errNevt[3]);
302     double corr_NevTTbar = hMC_[2][nvars_-1]->IntegralAndError(1,xbins[nvars_-1],errNevt[4]);
303     double corr_NevZJets = hMC_[3][nvars_-1]->IntegralAndError(1,xbins[nvars_-1],errNevt[5]);
304     double corr_NevSTtch = hMC_[4][nvars_-1]->IntegralAndError(1,xbins[nvars_-1],errNevt[6]);
305 jengbou 1.4 cout << "corr_NevSTtch = " << corr_NevSTtch << endl;
306     double corr_Nevmix = corr_NevQCD + corr_NevWJets + corr_NevTTbar
307     + corr_NevZJets + corr_NevSTtch;//should equal NevData for MC
308     errNevt[7] = sqrt(errNevt[1]*errNevt[1]+errNevt[3]*errNevt[3]+errNevt[4]*errNevt[4]+
309     errNevt[5]*errNevt[5]+errNevt[6]*errNevt[6]);
310 jengbou 1.3
311     // store nev in a vector
312     vNev_.push_back(NevData);
313     vNev_.push_back((useInv ? corr_NevQCD_NEW : corr_NevQCD));
314     vNev_.push_back(corr_NevWJets);
315    
316     // Non WJets (use MC expected values):
317     if (procTTbar > 0.) vNev_.push_back(corr_NevTTbar*procTTbar);
318     if (procZJets > 0.) vNev_.push_back(corr_NevZJets*procZJets);
319     if (procSTtch > 0.) vNev_.push_back(corr_NevSTtch*procSTtch);
320    
321 jengbou 1.4 outprint << "Events in Data = " << NevData << " +/- " << errNevt[0] << endl;
322     outprint << "Events QCD sample = " << corr_NevQCD << " +/- " << errNevt[1] << endl;
323     outprint << "Events InvSel sample = " << corr_NevQCD_NEW << " +/- " << errNevt[2] << endl;
324 jengbou 1.3 outprint << "---------------------------" << endl;
325 jengbou 1.4 outprint << "Events WJets sample = " << corr_NevWJets << " +/- " << errNevt[3] << endl;
326     outprint << "Events TTbar sample = " << corr_NevTTbar << " +/- " << errNevt[4] << endl;
327     outprint << "Events ZJets sample = " << corr_NevZJets << " +/- " << errNevt[5] << endl;
328     outprint << "Events STtch sample = " << corr_NevSTtch << " +/- " << errNevt[6] << endl;
329     outprint << "Total MC Events = " << corr_Nevmix << " +/- " << errNevt[7] << endl;
330 jengbou 1.3
331     //define the containers for chosen numbers (coressponding to the max KStest result)
332 jengbou 1.5 testMC maxProb[nvars_];
333 jengbou 1.3
334     //define the scale factors calculated using information obtained from all parameters
335     Double_t SFbackg = 0.0;
336     Double_t sumSFbackg = 0.0;
337     Double_t SFsample = 0.0;
338     Double_t sumSFsample = 0.0;
339     Double_t allKS = 0.0;
340    
341 jengbou 1.5 vector<testMC> resultsKS[nvars_];
342 jengbou 1.3 //do the KS test by varying the scale factors
343 jengbou 1.5 for (int i = 0; i < nvars_; i++) { // nvars_ variables
344 jengbou 1.3 TH1F *data = (TH1F*)hData_[i]->Clone("data");
345     data -> SetName("dataClone");
346     map<string,TH1F*> mHisto_;
347     mHisto_.clear();
348     mHisto_["Data"] = data;
349     mHisto_["QCD"] = (useInv ? (TH1F*)hQCD_NEW[i]->Clone() : (TH1F*)hMC_[0][i]->Clone());
350     mHisto_["WJets"] = (TH1F*)hMC_[1][i]->Clone();//WJets
351     if (procTTbar > 0.) {
352     mHisto_["TTbar"] = (TH1F*)hMC_[2][i]->Clone();
353     mHisto_["TTbar"]->Scale(procTTbar);
354     }
355     if (procZJets > 0.) {
356     mHisto_["ZJets"] = (TH1F*)hMC_[3][i]->Clone();
357     mHisto_["ZJets"]->Scale(procZJets);
358     }
359     if (procSTtch > 0.) {
360     mHisto_["STtch"] = (TH1F*)hMC_[4][i]->Clone();
361     mHisto_["STtch"]->Scale(procSTtch);
362     }
363 jengbou 1.2
364 jengbou 1.3 //data -> Scale(1./data->Integral());
365     resultsKS[i] = doKStest(vNev_,mHisto_);
366     testMC tksmax = getMax(resultsKS[i]);
367     maxProb[i] = tksmax;
368     outprint << "\nFor the plot " << histoLabelX[i] << " the results are:"<< endl;
369     outprint << "\tmax Probability = " << maxProb[i].prob << endl;
370     outprint << "\tproc_background = " << maxProb[i].scaleF_backg << endl;
371     outprint << "\tproc_sample = " << maxProb[i].scaleF_sample << endl;
372    
373     outprint << "\n\tpercent_B of Data = "
374     << maxProb[i].scaleF_backg*corr_NevQCD_NEW*100/NevData << endl;
375     outprint << "\tpercent_S of Data = "
376     << maxProb[i].scaleF_sample*corr_NevWJets*100/NevData << endl;
377     outprint << "---------------------------" << endl;
378    
379     //create the mixed samples with KS test results
380     sumSFbackg += maxProb[i].prob*maxProb[i].scaleF_backg;
381     sumSFsample += maxProb[i].prob*maxProb[i].scaleF_sample;
382     allKS += maxProb[i].prob;
383    
384     //fill a histogram with the results from the KS test for each variable
385     for (unsigned int jiter = 0; jiter < resultsKS[i].size(); jiter++) {
386     if (resultsKS[i].at(jiter).prob == 1.)
387     cout << "variable [" << i << "]: prob[" << jiter << "]= " << resultsKS[i].at(jiter).prob << endl;
388     hKSvalues_[i]->SetBinContent(jiter,resultsKS[i].at(jiter).prob);
389     }
390     delete data;
391 jengbou 1.2 }
392 jengbou 1.3 double SF_QCD = 0.;
393     //calculate the final scale factors
394     if (!useIKSopt) {
395     SFbackg = sumSFbackg/allKS;
396     SFsample = sumSFsample/allKS;
397     SF_QCD = SFbackg*corr_NevQCD_NEW/corr_NevQCD;
398     outprint << "allKS = " << allKS << "\tbackground = " << SFbackg << "\tsample = " << SFsample << endl;
399     outprint << "==> Scale Factor for QCD MC = " << SF_QCD << endl;
400     outprint << "\tcombined percent_B of Data = "
401     << SFbackg*corr_NevQCD_NEW*100/NevData << endl;
402     outprint << "\tcombined percent_S of Data = "
403     << SFsample*corr_NevWJets*100/NevData << endl;
404     outprint << "\n" << endl;
405     outprint << "=================================" << endl;
406     outprint << "\n" << endl;
407     } else {
408     testMC tksopt = getMaxL(resultsKS);
409     SFbackg = tksopt.scaleF_backg;
410     SFsample = tksopt.scaleF_sample;
411     SF_QCD = SFbackg*corr_NevQCD_NEW/corr_NevQCD;
412     outprint << "allKS (KSopt) = " << tksopt.prob
413     << "\tbackground = " << SFbackg << "\tsample = " << SFsample << endl;
414     outprint << "==> Scale Factor for QCD MC = " << SF_QCD << endl;
415     outprint << "\tcombined percent_B of Data = "
416     << SFbackg*corr_NevQCD_NEW*100/NevData << endl;
417     outprint << "\tcombined percent_S of Data = "
418     << SFsample*corr_NevWJets*100/NevData << endl;
419     outprint << "\n" << endl;
420     outprint << "=================================" << endl;
421     outprint << "\n" << endl;
422    
423     for (unsigned int jiter = 0; jiter < resultsKS[0].size(); jiter++) {
424 jengbou 1.5 double val_ = 1.;
425     for (size_t nvi = 0; nvi < nvars_; nvi++) {
426     val_ *= resultsKS[nvi].at(jiter).prob;
427     }
428     hKSvalues_[nvars_]->SetBinContent(jiter,val_);
429     //cout << "prod of IKS prob = " << val_ << endl;
430 jengbou 1.3 }
431     }
432 jengbou 1.1
433 jengbou 1.3 //=================================
434     // Plots
435     //=================================
436 jengbou 1.5 for (int i = 0; i < nvars_; i++) {// nvars_ variables
437 jengbou 1.3 gStyle->SetErrorX(0.);
438     hKSres_[i] -> Add((TH1F*)hQCD_NEW[i]->Clone(),(TH1F*)hMC_[1][i]->Clone(),SFbackg,SFsample);
439     hKSres_[i] -> Add((TH1F*)hMC_[2][i]->Clone(),procTTbar);
440     hKSres_[i] -> Add((TH1F*)hMC_[3][i]->Clone(),procZJets);
441     hKSres_[i] -> Add((TH1F*)hMC_[4][i]->Clone(),procSTtch);
442    
443     outprint << "hKSres->Integral() = " << hKSres_[i]->Integral() << endl;
444     outprint << "Data->Integral() = " << hData_[i]->Integral() << endl;
445    
446     hData_[i]->Rebin(2);
447     hQCD_NEW[i]->Rebin(2);
448     hMC_[0][i]->Rebin(2);
449     hMC_[1][i]->Rebin(2);
450     hKSres_[i]->Rebin(2);
451     //hKSvalues_[i]->Rebin(2);
452    
453     // Stack Wjets and QCD after scaled by KS factors
454     if (useInv && plotInvSelShape) {
455     hQCD_KS[i] = (TH1F*) hQCD_NEW[i]->Clone();
456     hQCD_KS[i]->Scale(SFbackg);
457     }
458     else {// Use MC QCD
459     hQCD_KS[i] = (TH1F*) hMC_[0][i]->Clone();
460     hQCD_KS[i]->Scale(SF_QCD);
461     }
462     hQCD_KS[i]->SetLineColor(style.QCDColor);
463     hQCD_KS[i]->SetFillColor(style.QCDColor);
464     hQCD_KS[i]->SetFillStyle(style.QCDFill);
465    
466     hWJets_KS[i] = (TH1F*) hMC_[1][i]->Clone();
467     hWJets_KS[i]->Scale(SFsample);
468     hWJets_KS[i]->SetLineColor(style.WJetsColor);
469     hWJets_KS[i]->SetFillColor(style.WJetsColor);
470     hWJets_KS[i]->SetFillStyle(style.WJetsFill);
471    
472     if (procTTbar > 0.) {
473     hMC_[2][i]->Rebin(2);
474     hMC_[2][i]->SetLineColor(style.TtbarColor);
475     hMC_[2][i]->SetFillColor(style.TtbarColor);
476     hMC_[2][i]->SetFillStyle(style.TtbarFill);
477     }
478     if (procZJets > 0.) {
479     hMC_[3][i]->Rebin(2);
480     hMC_[3][i]->SetLineColor(style.DYZJetsColor);
481     hMC_[3][i]->SetFillColor(style.DYZJetsColor);
482     hMC_[3][i]->SetFillStyle(style.DYZJetsFill);
483     }
484     if (procSTtch > 0.) {
485     hMC_[4][i]->Rebin(2);
486     hMC_[4][i]->SetLineColor(style.ST_tColor);
487     hMC_[4][i]->SetFillColor(style.ST_tColor);
488     hMC_[4][i]->SetFillStyle(style.ST_tFill);
489     }
490 jengbou 1.1
491 jengbou 1.3 THStack *hst = new THStack(invName,invName);
492     hst->Add((TH1F*)hQCD_KS[i]->Clone());
493     if (procSTtch > 0) {
494     TH1F *tmp = (TH1F*) hMC_[4][i]->Clone();
495     tmp->Scale(procSTtch);
496     hst->Add(tmp);
497     }
498     if (procZJets > 0) {
499     TH1F *tmp = (TH1F*) hMC_[3][i]->Clone();
500     tmp->Scale(procZJets);
501     hst->Add(tmp);
502     }
503     hst->Add((TH1F*)hWJets_KS[i]->Clone());
504     if (procTTbar > 0) {
505     TH1F *tmp = (TH1F*) hMC_[2][i]->Clone();
506     tmp->Scale(procTTbar);
507     hst->Add(tmp);
508     }
509 jengbou 1.1
510 jengbou 1.3 // Set plotting parameters
511     hData_[i] ->SetLineColor(1);
512     hQCD_NEW[i] ->SetLineColor(2);
513     hMC_[0][i] ->SetLineColor(4);
514     hMC_[1][i] ->SetLineColor(3);
515     hKSres_[i] ->SetLineColor(9);
516     hKSvalues_[i]->SetLineColor(i+1);
517    
518     hData_[i] ->SetMarkerColor(1);
519     hQCD_NEW[i] ->SetMarkerColor(2);
520     hMC_[0][i] ->SetMarkerColor(4);
521     hMC_[1][i] ->SetMarkerColor(3);
522     hKSres_[i] ->SetMarkerColor(9);
523     hKSvalues_[i]->SetMarkerColor(i+1);
524    
525     hData_[i] ->SetMarkerStyle((realData ? 20 : 24));
526     hQCD_NEW[i] ->SetMarkerStyle(20);
527     hMC_[0][i] ->SetMarkerStyle(20);
528     hMC_[1][i] ->SetMarkerStyle(20);
529     hKSres_[i] ->SetMarkerStyle(20);
530     hKSvalues_[i]->SetMarkerStyle(20);
531    
532     hData_[i] ->SetMarkerSize((realData ? 1.2 : 1.4));
533     hQCD_NEW[i] ->SetMarkerSize(1.1);
534     hMC_[0][i] ->SetMarkerSize(1.1);
535     hMC_[1][i] ->SetMarkerSize(1.1);
536     hKSres_[i] ->SetMarkerSize(1.18);
537     hKSvalues_[i]->SetMarkerSize(1.1);
538    
539     hData_[i] ->SetStats(0);
540     hQCD_NEW[i] ->SetStats(0);
541     hMC_[0][i] ->SetStats(0);
542     hMC_[1][i] ->SetStats(0);
543     hKSres_[i] ->SetStats(0);
544     hKSvalues_[i]->SetStats(0);
545     hQCD_KS[i] ->SetStats(0);
546     hWJets_KS[i] ->SetStats(0);
547    
548     hKSres_[i]->GetXaxis()->SetTitle(histoLabelX[i].c_str());
549     hKSres_[i]->GetYaxis()->SetTitle("Events");
550     hKSvalues_[i]->GetXaxis()->SetTitle("#lambda_{wjets}");
551     hKSvalues_[i]->GetYaxis()->SetTitle("KS test values");
552     hMC_[0][i]->GetXaxis()->SetTitle(histoLabelX[i].c_str());
553     hMC_[0][i]->GetYaxis()->SetTitle("A.U.");
554    
555     TLatex* text1 = new TLatex(1., 1., "2.88 pb^{-1} at #sqrt{s} = 7 TeV");
556     text1->SetNDC();
557     text1->SetTextAlign(13);
558     text1->SetX(0.63);
559     text1->SetY(0.91);
560     text1->SetTextFont(42);
561     text1->SetTextSize(0.03);
562     text1->SetTextSizePixels(20);// dflt=28
563    
564     // QCD shape comparison
565     TString cvsName1 = histoName[i]+"_QCD";
566     double val_Max_ = 0.;
567    
568     if(useInv) cvsName1 = cvsName1 + "_" + invName;
569 jengbou 1.4 cvsName1 += TString("_" + ssel_.str() + nthjbin_);
570 jengbou 1.5 cvs[cvsName1] = new TCanvas(cvsName1,"");
571 jengbou 1.3 hQCD_NEW[i] -> Scale(1./hQCD_NEW[i]->Integral());
572     hMC_[0][i] -> Scale(1./hMC_[0][i]->Integral());
573     hMC_[1][i] -> Scale(1./hMC_[1][i]->Integral());
574     outprint << "For " << histoName[i] << " , the KStest result btw MC_QCD/InvSel is = "
575     << hMC_[0][i] -> KolmogorovTest(hQCD_NEW[i],"") << endl;
576    
577     hMC_[0][i]->GetXaxis()->SetLabelFont(42);
578     hMC_[0][i]->GetXaxis()->SetLabelSize(0.04);
579     hMC_[0][i]->GetXaxis()->SetTitleFont(42);
580     hMC_[0][i]->GetXaxis()->SetTitleSize(0.05);
581     hMC_[0][i]->GetXaxis()->SetTitleOffset(1.2);
582     hMC_[0][i]->GetYaxis()->SetLabelFont(42);
583     hMC_[0][i]->GetYaxis()->SetLabelSize(0.04);
584     hMC_[0][i]->GetYaxis()->SetTitleFont(42);
585     hMC_[0][i]->GetYaxis()->SetTitleSize(0.05);
586     hMC_[0][i]->GetYaxis()->SetTitleOffset(1.3);
587     // Get better Y range
588     Int_t nbMax = hMC_[0][i]->GetMaximumBin();
589     double tmpval = hMC_[0][i]->GetBinContent(nbMax) + hMC_[0][i]->GetBinError(nbMax);
590     val_Max_ = (val_Max_ > 1.2*tmpval ? val_Max_ : 1.2*tmpval);
591 jengbou 1.1
592 jengbou 1.3 if (useInv) {
593     nbMax = hQCD_NEW[i]->GetMaximumBin();
594     tmpval = hQCD_NEW[i]->GetBinContent(nbMax) + hQCD_NEW[i]->GetBinError(nbMax);
595     val_Max_ = (val_Max_ > 1.2*tmpval ? val_Max_ : 1.2*tmpval);
596     }
597     hMC_[0][i]->GetYaxis()->SetRangeUser(0.,(val_Max_ > 1. ? 1. : val_Max_));
598     if (i == 0) // mu pT plot
599     hMC_[0][0]->GetXaxis()->SetRangeUser(20.,xhigh[0]);
600     hMC_[0][i]->DrawCopy("histe");
601    
602     if (useInv)
603     hQCD_NEW[i]->DrawCopy("samehiste");
604    
605     hMC_[1][i]->DrawCopy("samehiste");
606     TLegend *legend1 = new TLegend(0.6, 0.65, 0.92, 0.9);
607     legend1->SetFillColor(0);
608     legend1->SetFillStyle(0);
609    
610     legend1->AddEntry(hMC_[0][i], "QCD");
611     if (useInv)
612     if (!realData) legend1->AddEntry(hQCD_NEW[i], TString("QCD (" + invLabels[ninv] + ")"));
613     else legend1->AddEntry(hQCD_NEW[i], TString("Data (" + invLabels[ninv]) + ")");
614    
615     legend1->AddEntry(hMC_[1][i], style.WJetsText);
616     legend1->Draw();
617     legend1->SetFillColor(kWhite);
618     //latex->DrawLatex(0.22,0.91,histoName[i].c_str());
619     //cvs[cvsName1]->SetLogy();
620     TString nameCanvas1 = desDir+histoName[i]+"_QCD_"+suffix1_;
621     cvs[cvsName1]->SaveAs(nameCanvas1+".pdf");
622     cvs[cvsName1]->SaveAs(nameCanvas1+".png");
623    
624    
625     // KS results vs. (pseudo) Data
626     TString cvsName2 = histoName[i]+"_dataKS";
627     if(useInv) cvsName2 = cvsName2 + "_" + invName;
628 jengbou 1.4 cvsName2 += TString("_" + ssel_.str() + nthjbin_);
629 jengbou 1.5 cvs[cvsName2] = new TCanvas(cvsName2,"");
630 jengbou 1.3
631     // Get better Y range
632     val_Max_ = 0.;
633     nbMax = hData_[i]->GetMaximumBin();
634     tmpval = hData_[i]->GetBinContent(nbMax) + hData_[i]->GetBinError(nbMax);
635     val_Max_ = (val_Max_ > 1.05*tmpval ? val_Max_ : 1.05*tmpval);
636     nbMax = hKSres_[i]->GetMaximumBin();
637     tmpval = hKSres_[i]->GetBinContent(nbMax) + hKSres_[i]->GetBinError(nbMax);
638     val_Max_ = (val_Max_ > 1.05*tmpval ? val_Max_ : 1.05*tmpval);
639    
640     hData_[i]->GetYaxis()->SetRangeUser(0.,val_Max_);
641    
642     if (i == 0) // mu pT plot
643     hData_[0]->GetXaxis()->SetRangeUser(20.,xhigh[0]);
644     hData_[i]->Draw("P");
645     hst->Draw("samehist");
646     hKSres_[i]->Draw("samehiste");
647     hData_[i]->Draw("sameP");
648     TLegend *legend2 = new TLegend(0.7, 0.64, 0.9, 0.88);
649     legend2->AddEntry(hData_[i], (realData ? "Data" : "Pseudo Data"),"p");
650     legend2->AddEntry(hKSres_[i], "KS result");
651     if (procTTbar > 0.) legend2->AddEntry(hMC_[2][i], style.TtbarText);
652     legend2->AddEntry(hWJets_KS[i], style.WJetsText);
653     if (procZJets > 0.) legend2->AddEntry(hMC_[3][i], style.DYZJetsText);
654     if (procSTtch > 0.) legend2->AddEntry(hMC_[4][i], style.ST_tText);
655     if (realData)
656     legend2->AddEntry(hQCD_KS[i], ((!plotInvSelShape || !useInv) ? style.QCDText :
657     TString("#splitline{Data (" + invLabels[ninv] + ",}{Scaled by IKS)}")));
658     else
659     legend2->AddEntry(hQCD_KS[i], ((!plotInvSelShape || !useInv) ? style.QCDText :
660     TString("QCD (" + invLabels[ninv] + ")")));
661     legend2->SetFillColor(0);
662     legend2->SetFillStyle(0);
663     legend2->Draw();
664     legend2->SetFillColor(kWhite);
665     //latex->DrawLatex(0.22,0.91,histoName[i].c_str());
666     //cvs[cvsName2]->SetLogy();
667     text1->Draw();
668     TString nameCanvas2 = desDir+histoName[i]+"_dataKS_"+suffix1_;
669     cvs[cvsName2]->SaveAs(nameCanvas2+".pdf");
670     cvs[cvsName2]->SaveAs(nameCanvas2+".png");
671    
672     }
673    
674     // KS fit results
675     TString cvsName3 = "KStestValues";
676     if(useInv) cvsName3 = cvsName3 + "_" + invName;
677 jengbou 1.4 cvsName3 += TString("_" + ssel_.str() + nthjbin_);
678 jengbou 1.5 cvs[cvsName3] = new TCanvas(cvsName3,"");
679 jengbou 1.2 // Zoom-in
680 jengbou 1.3 Int_t nbins_ = hKSvalues_[0]->GetNbinsX();;
681     Int_t nbMax_ = 1;
682     Int_t nbMin_ = nbins_;
683     double val_Min_ = 1.;
684 jengbou 1.5 for (int i = 0; i < nvars_; ++i) {
685 jengbou 1.3 Int_t tmp = findFirstBinAbove(yLowBound, hKSvalues_[i]);
686     val_Min_ = (hKSvalues_[i]->GetBinContent(tmp) < 0.8*val_Min_ ?
687     hKSvalues_[i]->GetBinContent(tmp) : 0.8*val_Min_);
688     nbMin_ = (nbMin_ < tmp ? nbMin_ : tmp);
689     tmp = findLastBinAbove(0., hKSvalues_[i]);
690     nbMax_ = (nbMax_ > tmp ? nbMax_ : tmp);
691     }
692 jengbou 1.2 nbMin_ -= (0.1/stepsize);
693     nbMax_ += (0.1/stepsize);
694    
695 jengbou 1.3 hKSvalues_[0]->GetXaxis()->SetRange(nbMin_,nbMax_);
696     hKSvalues_[0]->GetYaxis()->SetRangeUser((1e-36 < val_Min_ ? val_Min_ : 1e-36),1.1);
697     hKSvalues_[0]->Draw();
698 jengbou 1.5 if (nvars_ > 1)
699     for (int i = 1; i < nvars_; i++) {
700     hKSvalues_[i]->Draw("same");
701     }
702    
703     TLegend *legend3 = new TLegend(0.2, (0.92-nvars_*0.07), 0.38, 0.92);
704     for (size_t i = 0; i < nvars_; i++)
705     legend3->AddEntry(hKSvalues_[i], TString(histoTitle[i]));
706    
707 jengbou 1.3 legend3->SetFillColor(0);
708     legend3->SetFillStyle(0);
709     legend3->Draw();
710     legend3->SetFillColor(kWhite);
711     //latex->DrawLatex(0.22,0.91,"KS test values");
712     cvs[cvsName3]->SetLogy();
713     TString nameCanvas3 = desDir+"KStestValues_newQCD_"+suffix1_;
714     cvs[cvsName3]->SaveAs(nameCanvas3+".pdf");
715     cvs[cvsName3]->SaveAs(nameCanvas3+".png");
716    
717     // KSopt
718     if (useIKSopt) {
719     TString cvsName4 = "KStestValues";
720     if(useInv) cvsName4 = cvsName4 + "_" + invName;
721 jengbou 1.4 cvsName4 += TString("_" + ssel_.str() + nthjbin_ + "_KSopt");
722 jengbou 1.3
723 jengbou 1.5 cvs[cvsName4] = new TCanvas(cvsName4,"");
724 jengbou 1.3 // Zoom-in
725 jengbou 1.5 nbins_ = hKSvalues_[nvars_]->GetNbinsX();;
726 jengbou 1.3 nbMax_ = 1;
727     nbMin_ = nbins_;
728     val_Min_ = 1.;
729 jengbou 1.5 Int_t tmp = findFirstBinAbove(1.e-50, hKSvalues_[nvars_]);
730     val_Min_ = (hKSvalues_[nvars_]->GetBinContent(tmp) < 0.8*val_Min_ ?
731     hKSvalues_[nvars_]->GetBinContent(tmp) : 0.8*val_Min_);
732 jengbou 1.3 nbMin_ = (nbMin_ < tmp ? nbMin_ : tmp);
733 jengbou 1.5 tmp = findLastBinAbove(0., hKSvalues_[nvars_]);
734 jengbou 1.3 nbMax_ = (nbMax_ > tmp ? nbMax_ : tmp);
735     nbMin_ -= (0.1/stepsize);
736     nbMax_ += (0.1/stepsize);
737    
738 jengbou 1.5 hKSvalues_[nvars_]->GetXaxis()->SetRange(nbMin_,nbMax_);
739     //hKSvalues_[nvars_]->GetYaxis()->SetRangeUser(1.e-100,1.e-50);
740     hKSvalues_[nvars_]->GetXaxis()->SetTitle("#lambda_{wjets}");
741     hKSvalues_[nvars_]->SetLineColor(4);
742     hKSvalues_[nvars_]->SetMarkerColor(4);
743     hKSvalues_[nvars_]->SetMarkerStyle(20);
744     hKSvalues_[nvars_]->SetMarkerSize(1.1);
745     hKSvalues_[nvars_]->SetStats(0);
746     hKSvalues_[nvars_]->Draw();
747 jengbou 1.3
748     TLegend *legend4 = new TLegend(0.2, 0.75, 0.6, 0.93);
749 jengbou 1.5 legend4->AddEntry(hKSvalues_[nvars_], "Product of KS test values");
750 jengbou 1.3 legend4->SetFillColor(0);
751     legend4->SetFillStyle(0);
752     legend4->Draw();
753     legend4->SetFillColor(kWhite);
754     cvs[cvsName4]->SetLogy();
755     TString nameCanvas4 = desDir+"KStestValues_newQCD_"+suffix1_+"_KSopt";
756     cvs[cvsName4]->SaveAs(nameCanvas4+".pdf");
757     cvs[cvsName4]->SaveAs(nameCanvas4+".png");
758     }
759     } // End of njbin loop
760     } // End of nsel loop
761 jengbou 1.1 } // End of ninv loop
762     }