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