ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/Jeng/scripts/ikstest.cc

Comparing UserCode/Jeng/scripts/ikstest.cc (file contents):
Revision 1.3 by jengbou, Thu Sep 2 05:36:21 2010 UTC vs.
Revision 1.8 by jengbou, Wed Sep 8 06:47:15 2010 UTC

#	Line 1 | Line 1
1	<	#include <vector>
2	<	#include <TFile.h>
3	<	#include <TH1.h>
4	<	#include <TH2.h>
5	<	#include <TLatex.h>
6	<	#include <TCanvas.h>
7	<	#include <TMath.h>
8	<	#include <TLegend.h>
9	<	#include <TObject.h>
10	<	#include <iostream>
11	<	#include <sstream>
12	<	#include <fstream>
13	<	#include <iomanip>
14	<	#include <map>
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	<	// User defined parameters
23	<	bool useInv           = false; // whether to use n-1 QCD template
24	<	bool realData         = false;
25	<	const double stepsize = 0.001;
26	<
27	<	//define the constants: 2.78/pb
24	<	const double weight_[3]    = {0.0506107, //QCD
25	<	0.0088222, //WJets
26	<	0.000295   //TTbar
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 procttjets  = 1.0;
32	<
33	<	struct testMC {
33	<	testMC(Double_t p = 0., Double_t sb = 0., Double_t ss = 0.){prob=p; scaleF_backg = sb; scaleF_sample = ss;}
34	<	Double_t prob;
35	<	Double_t scaleF_backg;
36	<	Double_t scaleF_sample;
37	<	};
38	<
39	<	//function for doing KS test
40	<	vector<testMC> doKStest(Double_t NsS, Double_t Ns1, Double_t Ns2, TH1F* mixS, TH1F* s1, TH1F* s2) {
41	<	vector<testMC> output;
42	<	//define the scale factors
43	<	Double_t sf1 = 0.0; // QCD
44	<	Double_t sf2 = 0.0; // Wjets
45	<	//KS test
46	<	do {
47	<	sf1 = (NsS - Ns2*sf2)/Ns1;
48	<	if (sf1 < 0) break;
49	<	//cout << "..........sf1 = " << sf1 << endl;
50	<	int nbins = mixS->GetNbinsX();
51	<	double xmin = mixS->GetXaxis()->GetBinLowEdge(1);
52	<	double xmax = mixS->GetXaxis()->GetBinUpEdge(nbins);
53	<	TH1F *test = new TH1F("test", "", nbins, xmin, xmax);
54	<	test -> Sumw2();
55	<	test -> Add(s1,s2,sf1,sf2);
56	<	//test->Scale(1./(1.*test->Integral()));
57	<	Double_t probability = mixS -> KolmogorovTest(test,"");
58	<	testMC temp = testMC(probability,sf1,sf2);
59	<	output.push_back(temp);
60	<	//     cout << "probability = " << setw(15) << temp.prob
61	<	//       << "; sfQCD = "     << setw(10) << temp.scaleF_backg
62	<	//       << "; sfWjets = "   << setw(6)  << temp.scaleF_sample << endl;
63	<	delete test;
64	<	sf2 = sf2 + stepsize;
65	<	} while(sf1 > 0 && sf2 <= 2.0);
66	<	return output;
67	<	}
68	<
69	<	//get the maximum KS test result
70	<	testMC getMax(vector<testMC> vec) {
71	<	testMC maxKSRes;
72	<	Double_t maximum = 0.0;
73	<	for (size_t i = 0; i < vec.size(); i++) {
74	<	if (maximum < vec.at(i).prob)  {
75	<	maximum = vec.at(i).prob;
76	<	maxKSRes = vec.at(i);
77	<	}
78	<	}
79	<	cout << "for maximum: " << setw(12) << maxKSRes.prob
80	<	<< "; sb = " << setw(10) << maxKSRes.scaleF_backg
81	<	<< "; ss = " << setw(5) << maxKSRes.scaleF_sample << endl;
82	<	return maxKSRes;
83	<	}
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	+	TString baseDir = "Results_2.88pb-1_NEW/";
37	+	// User defined parameters
38	+	bool useInv     = false; // whether to use n-1 QCD template
39	+	bool realData   = false;
40	+	// Ntuples to use
41	+	TString suffix  = "Sel0"; // Suffix of selection
42	+	TString invNames[2] = {"RelIsogt0p1","D0gt0p02"};
43	+	map<TString,TCanvas*> cvs; // map of usual histogram
44	+	bool debug_ = false;
45
46		//=================================
47		// Main program
48		//=================================
49		void ikstest() {
50	<	//Style();
50	>	CMSTopStyle style;
51	>	gROOT->SetStyle("CMS");
52		TLatex *latex = new TLatex();
53		latex->SetNDC();
54
55	<	ofstream outprint( "ikstest_results_20100901.txt" );
56	<	//open the files with histograms
57	<	map<string,TFile*> mfile;
58	<	mfile["Data"] = TFile::Open("skimmed_Data_20100901/Data_RefSel_v3.root");
59	<	// n-1 cuts
60	<	if (useInv) {
61	<	if (realData)
101	<	//       mfile["InvSel"] = TFile::Open("skimmed_Data_20100825/Data_D0ge0p03.root");
102	<	mfile["InvSel"] = TFile::Open("skimmed_Data_20100901/Data_RelIsoge0p1_v3.root");
103	<	else
104	<	//       mfile["InvSel"] = TFile::Open("skimmed_MC/QCD_D0ge0p03.root");
105	<	mfile["InvSel"] = TFile::Open("skimmed_MC/QCD_RelIsoge0p1_v3.root");
55	>	struct stat stbDir;
56	>	if (stat(baseDir,&stbDir)){
57	>	cout << "Creating folder : " << baseDir << endl;
58	>	if (mkdir(baseDir,0755) == -1){
59	>	std::cerr << "Error creating folder" << endl;
60	>	return;
61	>	}
62		}
63
64	<	// RefSel MC
65	<	mfile["0"] = TFile::Open("skimmed_MC/QCD_RefSel_v3.root");
66	<	mfile["1"] = TFile::Open("skimmed_MC/WJets_RefSel_v3.root");
67	<	mfile["2"] = TFile::Open("skimmed_MC/TTbar_RefSel_v3.root");
68	<
69	<	//define histograms and related parameters
70	<	string histoName[3] = {"h_mu_pt_calo","h_met_calo","h_mt_calo"};
71	<	string histoLabelX[3] = {"p_{T}^{good Muons}", "E_{T}^{#nu}", "m_{T}^{W}"};
72	<	Int_t xbins[3] = {20,20,40};
73	<	Double_t xlow[3] = {0.,0.,0.};
74	<	Double_t xhigh[3] = {100.,100.,200.};
119	<	string sample[3] = {"QCD","Wjets","ttjets"};
120	<
121	<	TH1F* h_[9];
122	<	TH1F* mixh_[3];
123	<	TH1F* hQCD_NEW[3];
124	<	TH1F* hKSres_[3];
125	<	TH1F* hKSvalues_[3];
126	<
127	<	//load the histograms from the root files
128	<	for (int i = 0; i < 3; i++) {// 3 variables
129	<	//cout << "file[" << i << "] : " << endl;
130	<	string nameNewHisto = "mix_"+histoName[i];
131	<	string nameNewHistoSFKS = "finalSF_"+histoName[i];
132	<	string nameNewHistoKSvalues = "KSvalues_"+histoLabelX[i];
133	<
134	<	mixh_[i] = new TH1F(nameNewHisto.c_str(),"",xbins[i],xlow[i],xhigh[i]);
135	<	hKSres_[i] = new TH1F(nameNewHistoSFKS.c_str(),"",xbins[i],xlow[i],xhigh[i]);
136	<	hKSvalues_[i] = new TH1F(nameNewHistoKSvalues.c_str(),"",2./stepsize, stepsize, 2.+stepsize);
137	<
138	<	if (!useInv) {//use QCD MC sample
139	<	hQCD_NEW[i] = (TH1F*) mfile["0"]->Get(TString(histoName[i]))->Clone();
140	<	hQCD_NEW[i] -> Scale(weight_[0]);
141	<	hQCD_NEW[i] -> SetName((histoName[i]).c_str());
64	>	int size_ninv = (useInv ? 2 : 1);
65	>	for (int ninv = 0;ninv < size_ninv; ++ninv) {
66	>	TString invName = invNames[ninv];
67	>	TString desDir;
68	>	// create output directory
69	>	if (!useInv) {
70	>	if (!realData) desDir = baseDir + "MC/";
71	>	else desDir = baseDir + "Data_MC/";
72	>	} else {
73	>	if (!realData) desDir = baseDir + "MC_"+invName+"/";
74	>	else desDir = baseDir + "Data_"+invName+"/";
75		}
76	<	else {
77	<	hQCD_NEW[i] = (TH1F*) mfile["InvSel"]->Get(TString(histoName[i]));
78	<	if (!realData) hQCD_NEW[i] -> Scale(weight_[0]);
79	<	hQCD_NEW[i] -> SetName((histoName[i]).c_str());
76	>	struct stat stDir;
77	>	if (!stat(desDir,&stDir)){
78	>	cout << "Output folder exists! Continue to overwrite it? (enter to continue; 'q' for quit)" << endl;
79	>	char incmd;
80	>	cin.get(incmd);
81	>	if (incmd == 'q') return;
82	>	} else {
83	>	cout << "Creating folder : " << desDir << endl;
84	>	if (mkdir(desDir,0755) == -1){
85	>	std::cerr << "Error creating folder" << endl;
86	>	return;
87	>	}
88		}
89
90	<	mixh_[i] -> Sumw2();
91	<	hKSres_[i] -> Sumw2();
92	<	hKSvalues_[i] -> Sumw2();
93	<	}
90	>	ofstream outprint(TString(desDir+"Results_"+suffix+".txt"));
91	>	//open the files with histograms
92	>	map<string,TFile*> mfile;
93	>	mfile.clear();
94	>	mfile["Data"] = TFile::Open(TString("skimmed_Data_2.88pb-1/Data_"+suffix+".root"));
95	>	// n-1 cuts
96	>	if (useInv) {
97	>	if (realData)
98	>	mfile["InvSel"] = TFile::Open(TString("skimmed_Data_2.88pb-1/Data_"+suffix+"_"+invName+".root"));
99	>	else
100	>	mfile["InvSel"] = TFile::Open(TString("skimmed_MC/v5/QCD_"+suffix+"_"+invName+".root"));
101	>	}
102	>	// RefSel MC
103	>	mfile["0"] = TFile::Open(TString("skimmed_MC/v5/QCD_"+suffix+".root"));
104	>	mfile["1"] = TFile::Open(TString("skimmed_MC/v5/WJets_"+suffix+".root"));
105	>	mfile["2"] = TFile::Open(TString("skimmed_MC/v5/TTbar_"+suffix+".root"));
106	>	mfile["3"] = TFile::Open(TString("skimmed_MC/v5/ZJets_"+suffix+".root"));
107	>	mfile["4"] = TFile::Open(TString("skimmed_MC/v5/STtch_"+suffix+".root"));
108	>
109	>	//define histograms and related parameters
110	>	string histoName[3] = {"h_mu_pt_calo","h_met_calo","h_mt_calo"};
111	>	string histoLabelX[3] = {"p_{T}^{#mu} [GeV/c]", "#slash{E}_{T} [GeV/c]", "M_{T}^{W} [GeV/c^{2}]"};
112	>	Int_t xbins[3] = {20,20,40};
113	>	Double_t xlow[3] = {0.,0.,0.};
114	>	Double_t xhigh[3] = {100.,100.,200.};
115	>	string sample[5] = {"QCD","WJets","TTbar","ZJets","STtch"};
116	>
117	>	TH1F* hMC_[5][3];   // MC histograms
118	>	TH1F* hData_[3];    // Data or mix MC
119	>	TH1F* hQCD_NEW[3];  // InvSel QCD shape
120	>	TH1F* hKSres_[3];
121	>	TH1F* hKSvalues_[3];
122	>	TH1F* hQCD_KS[3];
123	>	TH1F* hWJets_KS[3];
124	>
125	>	//load the histograms from the root files
126	>	for (int vi = 0; vi < 3; ++vi) {// 3 variables
127	>	//cout << "file[" << vi << "] : " << endl;
128	>	string nameNewHisto = "mix_"+histoName[vi];
129	>	string nameNewHistoSFKS = "finalSF_"+histoName[vi];
130	>	string nameNewHistoKSvalues = "KSvalues_"+histoLabelX[vi];
131	>
132	>	hData_[vi] = new TH1F(nameNewHisto.c_str(),"",xbins[vi],xlow[vi],xhigh[vi]);
133	>	hKSres_[vi] = new TH1F(nameNewHistoSFKS.c_str(),"",xbins[vi],xlow[vi],xhigh[vi]);
134	>	hKSvalues_[vi] = new TH1F(nameNewHistoKSvalues.c_str(),"",2./stepsize, stepsize, 2.+stepsize);
135	>
136	>	ostringstream ssc;
137	>	ssc << vi;
138	>
139	>	if (!useInv) {//use QCD MC sample
140	>	hQCD_NEW[vi] = (TH1F*) mfile["0"]->Get(TString(histoName[vi]))->Clone();
141	>	hQCD_NEW[vi] -> Scale(weight_[0]);
142	>	hQCD_NEW[vi] -> SetName(TString("InvSel_"+histoName[vi]+"_"+ssc.str()));
143	>	}
144	>	else {
145	>	hQCD_NEW[vi] = (TH1F*) mfile["InvSel"]->Get(TString(histoName[vi]))->Clone();
146	>	if (!realData) hQCD_NEW[vi] -> Scale(weight_[0]);
147	>	hQCD_NEW[vi] -> SetName(TString("InvSel_"+histoName[vi]+"_"+ssc.str()));
148	>	}
149	>	if (debug_) cout << "hQCD_NEW[" << vi << "] @ " << hQCD_NEW[vi] << endl;
150	>
151	>	hData_[vi] -> Sumw2();
152	>	hKSres_[vi] -> Sumw2();
153	>	hKSvalues_[vi] -> Sumw2();
154	>	}
155
156	<	for (int n = 0; n < 3; ++n) {// 3 MC samples
157	<	for (int ihisto = 0; ihisto < 3; ihisto++) {// 3 variables
158	<	//cout << "Variable[" << ihisto << "]" << endl;
159	<	string histo_name = histoName[ihisto]+sample[n];
160	<	ostringstream ss;
161	<	ss << n;
162	<	h_[n*3+ihisto] = (TH1F*) mfile[ss.str()]->Get(TString(histoName[ihisto]))->Clone();
163	<	h_[n*3+ihisto] -> Scale(weight_[n]);
164	<	h_[n*3+ihisto] -> SetName(histo_name.c_str());
156	>	for (int n = 0; n < 5; ++n) {// 3 MC samples
157	>	for (int ihisto = 0; ihisto < 3; ihisto++) {// 3 variables
158	>	//cout << "Variable[" << ihisto << "]" << endl;
159	>	string histo_name = histoName[ihisto]+sample[n];
160	>	ostringstream ss;
161	>	ss << n;
162	>	hMC_[n][ihisto] = (TH1F*) mfile[ss.str()]->Get(TString(histoName[ihisto]))->Clone();
163	>	if (debug_) {
164	>	cout << "File[" << n << "] @" << mfile[ss.str()]
165	>	<< "; histo[" << ihisto << "] @ " <<  mfile[ss.str()]->Get(TString(histoName[ihisto]))
166	>	<<"; hMC_[" << n << "][" << ihisto << "] raw evts = "
167	>	<< setw(12) << hMC_[n][ihisto]->Integral();
168	>	}
169	>	hMC_[n][ihisto] -> Scale(weight_[n]);
170	>	hMC_[n][ihisto] -> SetName(TString("MC_"+sample[n]+"_"+histoName[ihisto]));
171	>	if (debug_) cout << "; weighted num evts = " << setw(8) << hMC_[n][ihisto]->Integral() << endl;
172	>	}
173		}
164	–	}
174
175	<	//create the mixed samples = "data"
176	<	TCanvas *canvas0 = new TCanvas("Data","Data distributions");
177	<	canvas0->Divide(3,1);
178	<	for (int i = 0; i < 3; i++) {
179	<	canvas0->cd(i+1);
171	<	if (!realData) {
172	<	mixh_[i] -> Add(h_[i],h_[i+3], procQCD,procWjets);
173	<	//mixh_[i] -> Add(mixh_[i],h_[i+6], 1,procttjets);
174	<	//cout << "histo_name: " << mixh_[0]->GetNbinsX() << endl;
175	>	//create the mixed samples = "data"
176	>	TString cvsName0 = "Data";
177	>	if (useInv) {
178	>	cvsName0 += "_";
179	>	cvsName0 += invName;
180		}
181	<	else {
182	<	TH1F *htmp = (TH1F*) mfile["Data"]->Get(TString(histoName[i]));
183	<	mixh_[i] -> Add(htmp,1.);
181	>	cvs[cvsName0] = new TCanvas(cvsName0,"Data distributions",600,700);
182	>	cvs[cvsName0]->Divide(3,1);
183	>	for (int i = 0; i < 3; i++) {
184	>	cvs[cvsName0]->cd(i+1);
185	>	if (!realData) {
186	>	hData_[i] -> Add(hMC_[0][i],hMC_[1][i], procQCD,procWJets);
187	>	hData_[i] -> Add(hMC_[2][i], procTTbar);
188	>	hData_[i] -> Add(hMC_[3][i], procZJets);
189	>	hData_[i] -> Add(hMC_[4][i], procSTtch);
190	>	}
191	>	else {
192	>	TH1F *htmp = (TH1F*) mfile["Data"]->Get(TString(histoName[i]));
193	>	hData_[i] -> Add(htmp,1.);
194	>	}
195	>	hData_[i]->GetXaxis()->SetTitle(histoLabelX[i].c_str());
196	>	hData_[i]->GetYaxis()->SetTitle("Entries");
197	>	hData_[i]->DrawClone();
198		}
199	<	mixh_[i]->DrawClone();
200	<	}
201	<	canvas0->SaveAs("Data_distributions.pdf");
199	>	cvs[cvsName0]->SaveAs(TString(desDir+"Data_distributions.pdf"));
200	>
201	>	//Calculate num of events for each sample
202	>	vector<double> vNev_;
203	>
204	>	double NevData                = hData_[2]->Integral();
205	>	double corr_NevQCD            = hMC_[0][2]->Integral();
206	>	double corr_NevQCD_NEW        = hQCD_NEW[2]->Integral();
207	>	double corr_NevWJets          = hMC_[1][2]->Integral();
208	>	double corr_NevTTbar          = hMC_[2][2]->Integral();
209	>	double corr_NevZJets          = hMC_[3][2]->Integral();
210	>	double corr_NevSTtch          = hMC_[4][2]->Integral();
211	>	cout << "corr_NevSTtch = " << corr_NevSTtch << endl;
212	>	//     double corr_Nevmix            = procQCD*corr_NevQCD + procWJets*corr_NevWJets
213	>	//       + procTTbar*corr_NevTTbar+procZJets*corr_NevZJets + procSTtch*corr_NevSTtch;//should equal NevData for MC
214	>	// store nev in a vector
215	>	vNev_.push_back(NevData);
216	>	vNev_.push_back((useInv ? corr_NevQCD_NEW : corr_NevQCD));
217	>	vNev_.push_back(corr_NevWJets);
218	>
219	>	// Non WJets (use MC expected values):
220	>	if (procTTbar > 0.) vNev_.push_back(corr_NevTTbar*procTTbar);
221	>	if (procZJets > 0.) vNev_.push_back(corr_NevZJets*procZJets);
222	>	if (procSTtch > 0.) vNev_.push_back(corr_NevSTtch*procSTtch);
223
184	–	//define the weight corrections for each sample
185	–	double NevData                = mixh_[2]->Integral();
186	–	double corr_NevQCD            = h_[2]->Integral();
187	–	double corr_NevQCD_NEW        = hQCD_NEW[2]->Integral();
188	–	double corr_NevWjets          = h_[5]->Integral();
189	–	double corr_Nevttjets         = h_[8]->Integral();
190	–	double corr_Nevmix            = procQCD*corr_NevQCD+procWjets*corr_NevWjets;
191	–	//double corr_Nevmix            = procQCD*corr_NevQCD+procWjets*corr_NevWjets+procttjets*corr_Nevttjets;
192	–	if (!realData)
193	–	outprint << "Events mix sample    = " << corr_Nevmix << endl;
194	–	else
224		outprint << "Events in Data       = " << NevData << endl;
225	<	outprint << "Events QCD sample    = " << corr_NevQCD << endl;
226	<	outprint << "Events Wjets sample  = " << corr_NevWjets << endl;
198	<	outprint << "Events InvSel sample = " << corr_NevQCD_NEW << endl;
199	<
200	<	//define the containers for chosen numbers (coressponding to the max KStest result)
201	<	testMC maxProb[3];
202	<
203	<	//define the scale factors calculated using information obtained from all parameters
204	<	Double_t SFbackg = 0.0;
205	<	Double_t sumSFbackg = 0.0;
206	<	Double_t SFsample = 0.0;
207	<	Double_t sumSFsample = 0.0;
208	<	Double_t allKS = 0.0;
209	<
210	<	//do the KS test by varying the scale factors
211	<	for (int i = 0; i < 3; i++) { // 3 variables
212	<	TH1F *data = (TH1F*)mixh_[i]->Clone();
213	<	data -> SetName("dataClone");
214	<	//data -> Scale(1./data->Integral());
215	<	vector<testMC> resultsKS = doKStest((realData ? NevData : corr_Nevmix),
216	<	(useInv ? corr_NevQCD_NEW : corr_NevQCD),
217	<	corr_NevWjets,
218	<	data, hQCD_NEW[i], h_[i+3]);
219	<	testMC tksmax = getMax(resultsKS);
220	<	maxProb[i] = tksmax;
221	<	outprint << "\nFor the plot " << histoLabelX[i] << " the results are:"<< endl;
222	<	outprint << "\tmax Probability = " << maxProb[i].prob << endl;
223	<	outprint << "\tproc_background = " << maxProb[i].scaleF_backg << endl;
224	<	outprint << "\tproc_sample     = " << maxProb[i].scaleF_sample << endl;
225	<
226	<	outprint << "\n\tpercent_B of Data = "
227	<	<< maxProb[i].scaleF_backg*corr_NevQCD_NEW*100/(realData ? NevData : corr_Nevmix) << endl;
228	<	outprint << "\tpercent_S of Data = "
229	<	<< maxProb[i].scaleF_sample*corr_NevWjets*100/(realData ? NevData : corr_Nevmix) << endl;
225	>	outprint << "Events QCD sample    = " << corr_NevQCD << endl;
226	>	outprint << "Events InvSel sample = " << corr_NevQCD_NEW << endl;
227		outprint << "---------------------------" << endl;
228	<
229	<	//create the mixed samples with KS test results
230	<	sumSFbackg += maxProb[i].prob*maxProb[i].scaleF_backg;
231	<	sumSFsample += maxProb[i].prob*maxProb[i].scaleF_sample;
232	<	allKS += maxProb[i].prob;
233	<
234	<	//fill a histogram with the results from the KS test for each variable
235	<	for (int jiter = 0; jiter < resultsKS.size(); jiter++) {
236	<	if (resultsKS.at(jiter).prob == 1.)
237	<	cout << "variable [" << i << "]: prob[" << jiter << "]= " << resultsKS.at(jiter).prob << endl;
238	<	hKSvalues_[i]->SetBinContent(jiter,resultsKS.at(jiter).prob);
228	>	outprint << "Events WJets sample  = " << corr_NevWJets << endl;
229	>	outprint << "Events TTbar sample  = " << corr_NevTTbar << endl;
230	>	outprint << "Events ZJets sample  = " << corr_NevZJets << endl;
231	>	outprint << "Events STtch sample  = " << corr_NevSTtch << endl;
232	>
233	>	//define the containers for chosen numbers (coressponding to the max KStest result)
234	>	testMC maxProb[3];
235	>
236	>	//define the scale factors calculated using information obtained from all parameters
237	>	Double_t SFbackg = 0.0;
238	>	Double_t sumSFbackg = 0.0;
239	>	Double_t SFsample = 0.0;
240	>	Double_t sumSFsample = 0.0;
241	>	Double_t allKS = 0.0;
242	>
243	>	//do the KS test by varying the scale factors
244	>	for (int i = 0; i < 3; i++) { // 3 variables
245	>	TH1F *data = (TH1F*)hData_[i]->Clone("data");
246	>	data -> SetName("dataClone");
247	>	map<string,TH1F*> mHisto_;
248	>	mHisto_.clear();
249	>	mHisto_["Data"]  = data;
250	>	mHisto_["QCD"]   = (useInv ? (TH1F*)hQCD_NEW[i]->Clone() : (TH1F*)hMC_[0][i]->Clone());
251	>	mHisto_["WJets"] = (TH1F*)hMC_[1][i]->Clone();//WJets
252	>	if (procTTbar > 0.) mHisto_["TTbar"] = (TH1F*)hMC_[2][i]->Clone();//TTbar
253	>	if (procZJets > 0.) mHisto_["ZJets"] = (TH1F*)hMC_[3][i]->Clone();//ZJets
254	>	if (procSTtch > 0.) mHisto_["STtch"] = (TH1F*)hMC_[4][i]->Clone();//STtch
255	>
256	>	//data -> Scale(1./data->Integral());
257	>	vector<testMC> resultsKS = doKStest(vNev_,mHisto_);
258	>	testMC tksmax = getMax(resultsKS);
259	>	maxProb[i] = tksmax;
260	>	outprint << "\nFor the plot " << histoLabelX[i] << " the results are:"<< endl;
261	>	outprint << "\tmax Probability = " << maxProb[i].prob << endl;
262	>	outprint << "\tproc_background = " << maxProb[i].scaleF_backg << endl;
263	>	outprint << "\tproc_sample     = " << maxProb[i].scaleF_sample << endl;
264	>
265	>	outprint << "\n\tpercent_B of Data = "
266	>	<< maxProb[i].scaleF_backg*corr_NevQCD_NEW*100/NevData << endl;
267	>	outprint << "\tpercent_S of Data = "
268	>	<< maxProb[i].scaleF_sample*corr_NevWJets*100/NevData << endl;
269	>	outprint << "---------------------------" << endl;
270	>
271	>	//create the mixed samples with KS test results
272	>	sumSFbackg += maxProb[i].prob*maxProb[i].scaleF_backg;
273	>	sumSFsample += maxProb[i].prob*maxProb[i].scaleF_sample;
274	>	allKS += maxProb[i].prob;
275	>
276	>	//fill a histogram with the results from the KS test for each variable
277	>	for (unsigned int jiter = 0; jiter < resultsKS.size(); jiter++) {
278	>	if (resultsKS.at(jiter).prob == 1.)
279	>	cout << "variable [" << i << "]: prob[" << jiter << "]= " << resultsKS.at(jiter).prob << endl;
280	>	hKSvalues_[i]->SetBinContent(jiter,resultsKS.at(jiter).prob);
281	>	}
282	>	delete data;
283		}
243	–	delete data;
244	–	}
284
285	<	//calculate the final scale factors
286	<	SFbackg = sumSFbackg/allKS;
287	<	SFsample = sumSFsample/allKS;
288	<	outprint << "allKS = " << allKS << "\tbackground = " << SFbackg << "\tsample = " << SFsample << endl;
289	<	outprint << "==> Scale Factor for QCD MC = " << SFbackg*corr_NevQCD_NEW/corr_NevQCD << endl;
290	<	outprint << "\tcombined percent_B of Data = "
291	<	<< SFbackg*corr_NevQCD_NEW*100/(realData ? NevData : corr_Nevmix) << endl;
292	<	outprint << "\tcombined percent_S of Data = "
293	<	<< SFsample*corr_NevWjets*100/(realData ? NevData : corr_Nevmix) << endl;
294	<	outprint << "\n" << endl;
295	<	outprint << "=================================" << endl;
296	<	outprint << "\n" << endl;
297	<
298	<
299	<	//=================================
300	<	// Plots
301	<	//=================================
302	<	for (int i = 0; i < 3; i++) {// 3 variables
303	<	hKSres_[i] -> Add(hQCD_NEW[i],h_[i+3],SFbackg,SFsample);
304	<	outprint << "hKSres->Integral() = " << hKSres_[i]->Integral() << endl;
305	<	outprint << "Data->Integral()   = " << mixh_[i]->Integral() << endl;
306	<
307	<	mixh_[i]->Rebin(2);
308	<	hQCD_NEW[i]->Rebin(2);
309	<	h_[i]->Rebin(2);
310	<	h_[i+3]->Rebin(2);
311	<	hKSres_[i]->Rebin(2);
312	<	//hKSvalues_[i]->Rebin(2);
313	<
314	<	mixh_[i]     ->SetLineColor(1);
315	<	hQCD_NEW[i]  ->SetLineColor(3);
316	<	h_[i]        ->SetLineColor(6);
317	<	h_[i+3]      ->SetLineColor(4);
318	<	hKSres_[i]   ->SetLineColor(2);
319	<	hKSvalues_[i]->SetLineColor(i+1);
320	<
321	<	mixh_[i]     ->SetMarkerColor(1);
322	<	hQCD_NEW[i]  ->SetMarkerColor(3);
323	<	h_[i]        ->SetMarkerColor(6);
324	<	h_[i+3]      ->SetMarkerColor(4);
325	<	hKSres_[i]   ->SetMarkerColor(2);
326	<	hKSvalues_[i]->SetMarkerColor(i+1);
327	<
328	<	mixh_[i]     ->SetMarkerStyle(24);
329	<	hQCD_NEW[i]  ->SetMarkerStyle(20);
330	<	h_[i]        ->SetMarkerStyle(20);
331	<	h_[i+3]      ->SetMarkerStyle(20);
332	<	hKSres_[i]   ->SetMarkerStyle(20);
333	<	hKSvalues_[i]->SetMarkerStyle(20);
334	<
335	<	mixh_[i]     ->SetMarkerSize(1.4);
336	<	hQCD_NEW[i]  ->SetMarkerSize(1.1);
337	<	h_[i]        ->SetMarkerSize(1.1);
338	<	h_[i+3]      ->SetMarkerSize(1.1);
339	<	hKSres_[i]   ->SetMarkerSize(0.9);
340	<	hKSvalues_[i]->SetMarkerSize(1.1);
341	<
342	<	mixh_[i]     ->SetStats(0);
343	<	hQCD_NEW[i]  ->SetStats(0);
344	<	h_[i]        ->SetStats(0);
345	<	h_[i+3]      ->SetStats(0);
346	<	hKSres_[i]   ->SetStats(0);
347	<	hKSvalues_[i]->SetStats(0);
348	<
349	<	mixh_[i]->GetXaxis()->SetTitle(histoLabelX[i].c_str());
350	<	mixh_[i]->GetYaxis()->SetTitle("Entries");
351	<	hKSres_[i]->GetXaxis()->SetTitle(histoLabelX[i].c_str());
352	<	hKSres_[i]->GetYaxis()->SetTitle("Entries");
353	<	hKSvalues_[i]->GetXaxis()->SetTitle("iteration #");
354	<	hKSvalues_[i]->GetYaxis()->SetTitle("KS test values");
355	<
356	<	string nameCanvas1 = histoName[i]+"_QCD.pdf";
357	<	TCanvas *canvas1 = new TCanvas(nameCanvas1.c_str(), "");
358	<	hQCD_NEW[i] -> Scale(1./hQCD_NEW[i]->Integral());
359	<	h_[i] -> Scale(1./h_[i]->Integral());
360	<	h_[i+3] -> Scale(1./h_[i+3]->Integral());
361	<	outprint << "For " << histoName[i] << " , the KStest result btw MC_QCD/InvSel is = "
362	<	<< h_[i] -> KolmogorovTest(hQCD_NEW[i],"") << endl;
363	<	hQCD_NEW[i]->Draw("P");
364	<	h_[i]->Draw("sameP");
365	<	h_[i+3]->Draw("sameP");
366	<	TLegend *legend1 = new TLegend(0.7, 0.70, 0.9, 0.85);
367	<	legend1->AddEntry(h_[i], "default");
368	<	legend1->AddEntry(h_[i+3], "W+jets");
369	<	legend1->AddEntry(hQCD_NEW[i], "new");
370	<	legend1->Draw();
371	<	legend1->SetFillColor(kWhite);
372	<	latex->DrawLatex(0.22,0.91,histoName[i].c_str());
373	<	canvas1->SetLogy();
374	<	canvas1->SaveAs(nameCanvas1.c_str());
375	<
376	<	string nameCanvas2 = histoName[i]+"_dataKS.pdf";
377	<	TCanvas *canvas2 = new TCanvas(nameCanvas2.c_str(), "");
378	<	hKSres_[i]->Draw("P");
379	<	mixh_[i]->Draw("sameP");
380	<	TLegend *legend2 = new TLegend(0.7, 0.70, 0.9, 0.85);
381	<	legend2->AddEntry(mixh_[i], "Data");
382	<	legend2->AddEntry(hKSres_[i], "KS result");
383	<	legend2->Draw();
384	<	legend2->SetFillColor(kWhite);
385	<	latex->DrawLatex(0.22,0.91,histoName[i].c_str());
386	<	canvas2->SetLogy();
387	<	canvas2->SaveAs(nameCanvas2.c_str());
388	<
389	<	}
285	>	//calculate the final scale factors
286	>	SFbackg = sumSFbackg/allKS;
287	>	SFsample = sumSFsample/allKS;
288	>	outprint << "allKS = " << allKS << "\tbackground = " << SFbackg << "\tsample = " << SFsample << endl;
289	>	outprint << "==> Scale Factor for QCD MC = " << SFbackg*corr_NevQCD_NEW/corr_NevQCD << endl;
290	>	outprint << "\tcombined percent_B of Data = "
291	>	<< SFbackg*corr_NevQCD_NEW*100/NevData << endl;
292	>	outprint << "\tcombined percent_S of Data = "
293	>	<< SFsample*corr_NevWJets*100/NevData << endl;
294	>	outprint << "\n" << endl;
295	>	outprint << "=================================" << endl;
296	>	outprint << "\n" << endl;
297	>
298	>
299	>	//=================================
300	>	// Plots
301	>	//=================================
302	>	for (int i = 0; i < 3; i++) {// 3 variables
303	>	hKSres_[i] -> Add((TH1F*)hQCD_NEW[i]->Clone(),(TH1F*)hMC_[1][i]->Clone(),SFbackg,SFsample);
304	>	hKSres_[i] -> Add((TH1F*)hMC_[2][i]->Clone(),procTTbar);
305	>	hKSres_[i] -> Add((TH1F*)hMC_[3][i]->Clone(),procZJets);
306	>	hKSres_[i] -> Add((TH1F*)hMC_[4][i]->Clone(),procSTtch);
307	>
308	>	outprint << "hKSres->Integral() = " << hKSres_[i]->Integral() << endl;
309	>	outprint << "Data->Integral()   = " << hData_[i]->Integral() << endl;
310	>
311	>	hData_[i]->Rebin(2);
312	>	hQCD_NEW[i]->Rebin(2);
313	>	hMC_[0][i]->Rebin(2);
314	>	hMC_[1][i]->Rebin(2);
315	>	hKSres_[i]->Rebin(2);
316	>	//hKSvalues_[i]->Rebin(2);
317	>
318	>	// Stack Wjets and QCD after scaled by KS factors
319	>	hQCD_KS[i] = (TH1F*) hQCD_NEW[i]->Clone();
320	>	hQCD_KS[i]->Scale(SFbackg);
321	>	hQCD_KS[i]->SetLineColor(style.QCDColor);
322	>	hQCD_KS[i]->SetFillColor(style.QCDColor);
323	>	hQCD_KS[i]->SetFillStyle(style.QCDFill);
324	>
325	>	hWJets_KS[i] = (TH1F*) hMC_[1][i]->Clone();
326	>	hWJets_KS[i]->Scale(SFsample);
327	>	hWJets_KS[i]->SetLineColor(style.WJetsColor);
328	>	hWJets_KS[i]->SetFillColor(style.WJetsColor);
329	>	hWJets_KS[i]->SetFillStyle(style.WJetsFill);
330	>
331	>	if (procTTbar > 0.) {
332	>	hMC_[2][i]->Rebin(2);
333	>	hMC_[2][i]->SetLineColor(style.TtbarColor);
334	>	hMC_[2][i]->SetFillColor(style.TtbarColor);
335	>	hMC_[2][i]->SetFillStyle(style.TtbarFill);
336	>	}
337	>	if (procZJets > 0.) {
338	>	hMC_[3][i]->Rebin(2);
339	>	hMC_[3][i]->SetLineColor(style.DYZJetsColor);
340	>	hMC_[3][i]->SetFillColor(style.DYZJetsColor);
341	>	hMC_[3][i]->SetFillStyle(style.DYZJetsFill);
342	>	}
343	>	if (procSTtch > 0.) {
344	>	hMC_[4][i]->Rebin(2);
345	>	hMC_[4][i]->SetLineColor(style.ST_t_sColor);
346	>	hMC_[4][i]->SetFillColor(style.ST_t_sColor);
347	>	hMC_[4][i]->SetFillStyle(style.ST_t_sFill);
348	>	}
349	>	THStack *hst = new THStack(invName,invName);
350	>	hst->Add((TH1F*)hQCD_KS[i]->Clone());
351	>	if (procSTtch > 0) hst->Add((TH1F*)hMC_[4][i]->Clone());
352	>	if (procZJets > 0) hst->Add((TH1F*)hMC_[3][i]->Clone());
353	>	hst->Add((TH1F*)hWJets_KS[i]->Clone());
354	>	if (procTTbar > 0) hst->Add((TH1F*)hMC_[2][i]->Clone());
355	>
356	>	// Set plotting parameters
357	>	hData_[i]    ->SetLineColor(1);
358	>	hQCD_NEW[i]  ->SetLineColor(2);
359	>	hMC_[0][i]   ->SetLineColor(4);
360	>	hMC_[1][i]   ->SetLineColor(3);
361	>	hKSres_[i]   ->SetLineColor(2);
362	>	hKSvalues_[i]->SetLineColor(i+1);
363	>
364	>	hData_[i]    ->SetMarkerColor(1);
365	>	hQCD_NEW[i]  ->SetMarkerColor(2);
366	>	hMC_[0][i]   ->SetMarkerColor(4);
367	>	hMC_[1][i]   ->SetMarkerColor(3);
368	>	hKSres_[i]   ->SetMarkerColor(2);
369	>	hKSvalues_[i]->SetMarkerColor(i+1);
370	>
371	>	hData_[i]    ->SetMarkerStyle(24);
372	>	hQCD_NEW[i]  ->SetMarkerStyle(20);
373	>	hMC_[0][i]   ->SetMarkerStyle(20);
374	>	hMC_[1][i]   ->SetMarkerStyle(20);
375	>	hKSres_[i]   ->SetMarkerStyle(20);
376	>	hKSvalues_[i]->SetMarkerStyle(20);
377	>
378	>	hData_[i]    ->SetMarkerSize(1.4);
379	>	hQCD_NEW[i]  ->SetMarkerSize(1.1);
380	>	hMC_[0][i]   ->SetMarkerSize(1.1);
381	>	hMC_[1][i]   ->SetMarkerSize(1.1);
382	>	hKSres_[i]   ->SetMarkerSize(0.9);
383	>	hKSvalues_[i]->SetMarkerSize(1.1);
384	>
385	>	hData_[i]    ->SetStats(0);
386	>	hQCD_NEW[i]  ->SetStats(0);
387	>	hMC_[0][i]   ->SetStats(0);
388	>	hMC_[1][i]   ->SetStats(0);
389	>	hKSres_[i]   ->SetStats(0);
390	>	hKSvalues_[i]->SetStats(0);
391	>	hQCD_KS[i]   ->SetStats(0);
392	>	hWJets_KS[i] ->SetStats(0);
393	>
394	>	hKSres_[i]->GetXaxis()->SetTitle(histoLabelX[i].c_str());
395	>	hKSres_[i]->GetYaxis()->SetTitle("Entries");
396	>	hKSvalues_[i]->GetXaxis()->SetTitle("iteration #");
397	>	hKSvalues_[i]->GetYaxis()->SetTitle("KS test values");
398	>	hMC_[0][i]->GetXaxis()->SetTitle(histoLabelX[i].c_str());
399	>	hMC_[0][i]->GetYaxis()->SetTitle("A.U.");
400	>
401	>
402	>	TString nameCanvas1 = desDir+histoName[i]+"_QCD_"+suffix+".pdf";
403	>	TString cvsName1 = histoName[i]+"_QCD";
404	>	if(useInv) cvsName1 = cvsName1 + "_" + invName;
405	>	cvs[cvsName1] = new TCanvas(cvsName1,"",600,700);
406	>	hQCD_NEW[i] -> Scale(1./hQCD_NEW[i]->Integral());
407	>	hMC_[0][i] -> Scale(1./hMC_[0][i]->Integral());
408	>	hMC_[1][i] -> Scale(1./hMC_[1][i]->Integral());
409	>	outprint << "For " << histoName[i] << " , the KStest result btw MC_QCD/InvSel is = "
410	>	<< hMC_[0][i] -> KolmogorovTest(hQCD_NEW[i],"") << endl;
411	>	hMC_[0][i]->DrawCopy("P");
412	>	if (useInv)
413	>	hQCD_NEW[i]->DrawCopy("sameP");
414	>	hMC_[1][i]->DrawCopy("sameP");
415	>	TLegend *legend1 = new TLegend(0.7, 0.65, 0.9, 0.85);
416	>	legend1->AddEntry(hMC_[0][i], "QCD");
417	>	if (useInv)
418	>	legend1->AddEntry(hQCD_NEW[i], "QCD - InvSel");
419	>	legend1->AddEntry(hMC_[1][i], style.WJetsText);
420	>	legend1->Draw();
421	>	legend1->SetFillColor(kWhite);
422	>	//latex->DrawLatex(0.22,0.91,histoName[i].c_str());
423	>	//cvs[cvsName1]->SetLogy();
424	>	cvs[cvsName1]->SaveAs(nameCanvas1);
425	>
426	>
427	>	TString nameCanvas2 = desDir+histoName[i]+"_dataKS_"+suffix+".pdf";
428	>	TString cvsName2 = histoName[i]+"_dataKS";
429	>	if(useInv) cvsName2 = cvsName2 + "_" + invName;
430	>	cvs[cvsName2] = new TCanvas(cvsName2,"",600,700);
431	>	hst->Draw("hist");
432	>	hKSres_[i]->Draw("sameP");
433	>	hData_[i]->Draw("sameP");
434	>
435	>	TLegend *legend2 = new TLegend(0.7, 0.65, 0.9, 0.85);
436	>	legend2->AddEntry(hData_[i], "Data");
437	>	legend2->AddEntry(hKSres_[i], "KS result");
438	>	if (procTTbar > 0.) legend2->AddEntry(hMC_[2][i], style.TtbarText);
439	>	legend2->AddEntry(hWJets_KS[i], style.WJetsText);
440	>	if (procZJets > 0.) legend2->AddEntry(hMC_[3][i], style.DYZJetsText);
441	>	if (procSTtch > 0.) legend2->AddEntry(hMC_[4][i], style.ST_t_sText);
442	>	legend2->AddEntry(hQCD_KS[i], style.QCDText);
443	>	legend2->Draw();
444	>	legend2->SetFillColor(kWhite);
445	>	//latex->DrawLatex(0.22,0.91,histoName[i].c_str());
446	>	//cvs[cvsName2]->SetLogy();
447	>	cvs[cvsName2]->SaveAs(nameCanvas2);
448
449	+	}
450
451	<	TCanvas *canvas3 = new TCanvas("KStestValues", "");
452	<	//hKSvalues_[0]->GetXaxis()->SetRangeUser(0.9,1.1);
453	<	//hKSvalues_[0]->GetYaxis()->SetRangeUser(1e-3,1.1);
454	<	hKSvalues_[0]->Draw();
455	<	hKSvalues_[1]->SetLineColor(2);
456	<	hKSvalues_[1]->Draw("same");
457	<	hKSvalues_[2]->SetLineColor(4);
458	<	hKSvalues_[2]->Draw("same");
459	<	TLegend *legend3 = new TLegend(0.7, 0.70, 0.9, 0.85);
460	<	legend3->AddEntry(hKSvalues_[0], "muon_pT");
461	<	legend3->AddEntry(hKSvalues_[1], "MET");
462	<	legend3->AddEntry(hKSvalues_[2], "W_mT");
463	<	legend3->Draw();
464	<	legend3->SetFillColor(kWhite);
465	<	latex->DrawLatex(0.22,0.91,"KS test values");
466	<	canvas3->SetLogy();
467	<	string nameCanvas3 = "KStestValues_newQCD.pdf";
468	<	canvas3->SaveAs(nameCanvas3.c_str());
469	<
451	>	TString cvsName3 = "KStestValues";
452	>	if(useInv) cvsName3 = cvsName3 + "_" + invName;
453	>	cvs[cvsName3] = new TCanvas(cvsName3,"",600,700);
454	>	if (!realData) hKSvalues_[0]->GetXaxis()->SetRangeUser(0.95,1.1);
455	>	hKSvalues_[0]->GetYaxis()->SetRangeUser(1e-36,1.2);
456	>	hKSvalues_[0]->Draw();
457	>	hKSvalues_[1]->Draw("same");
458	>	hKSvalues_[2]->Draw("same");
459	>	TLegend *legend3 = new TLegend(0.7, 0.70, 0.9, 0.85);
460	>	legend3->AddEntry(hKSvalues_[0], "muon_pT");
461	>	legend3->AddEntry(hKSvalues_[1], "MET");
462	>	legend3->AddEntry(hKSvalues_[2], "W_mT");
463	>	legend3->Draw();
464	>	legend3->SetFillColor(kWhite);
465	>	//latex->DrawLatex(0.22,0.91,"KS test values");
466	>	cvs[cvsName3]->SetLogy();
467	>	TString nameCanvas3 = desDir+"KStestValues_newQCD"+suffix+".pdf";
468	>	cvs[cvsName3]->SaveAs(nameCanvas3);
469	>	}
470		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines