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

File Contents

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