ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/Jeng/scripts/fitRelIso_Data.C

Comparing UserCode/Jeng/scripts/fitRelIso_Data.C (file contents):
Revision 1.1 by jengbou, Thu Jun 24 02:16:25 2010 UTC vs.
Revision 1.3 by jengbou, Wed Jul 21 00:55:58 2010 UTC

#	Line 19 | Line 19
19		#include "map.h"
20		#include <sstream>
21		#include <fstream>
22	+	#include "TopStyle/tdrstyle.C"
23	+	#include "TopStyle/CMSTopStyle.cc"
24
25		using namespace std;
26
#	Line 27 | Line 29 | double sf_mc = 1.;
29		// Dynamic(true) or static(false) range
30		bool dynamic      = true;
31		TString isoname   = "New";
32	<	TString s0        = "fuser";
33	<	TString fitOpt    = "QLM0";
32	>	TString s0        = "fuser";//fuser,Landau,gauss
33	>	TString fitOpt    = "QLLMER0";
34		bool debug        = true;
35		double yMax       = 0.;
36		bool DrawOverFlow = true;
37	<	bool drawtail     = false;
37	>	bool drawtail     = true;
38		double wxmin      = 1.;
39		double wxmax      = 1.;
40		double plot_xMax  = 1.4; // max of x axis for output plots; 1.0 for Old reliso
41	+	double plot_xMin  = 0.3; // min of x axis for for Old reliso
42	+	int scaleType     = 0; // 0: scale MC w.r.t. total number of events; 1: scale QCD number of event to fitted result
43	+	bool drawStats    = false;
44
45		map<TString,TCanvas*> cvs; // map of usual histogram
46
#	Line 78 | Line 83 | void cmsPrel(double intLumi){
83
84		void fitRelIso_Data()
85		{
86	+	//setTDRStyle();
87	+	CMSTopStyle style;
88	+	gROOT->SetStyle("CMS");
89	+	style.setupICHEPv1();
90	+
91		//New RelIso
92		if (isoname == "New") {
93	<	s0 = "landau";
94	<	wxmin = 10.;
95	<	wxmax = wxmin/3.;
93	>	//s0 = "landau";
94	>	wxmax = 4.;//4
95	>	wxmin = 2.;
96	>	plot_xMin = 0.;
97		// Select fitting function:
98		}
99		//Old RelIso
100		if (isoname == "Old") {
101		s0 = "gaus"; // for old
102	<	wxmin = 1.;
103	<	wxmax = 1./3.;
102	>	//s0 = "fuser";
103	>	wxmin = 2.;
104	>	wxmax = 1.;
105		DrawOverFlow = false;
106		plot_xMax = 1.;
107		}
108
97	–	gROOT->SetStyle("CMS");
98	–
109		cout << " =========================== " << endl;
110		cout << " = Start fitting           = " << endl;
111		cout << " =========================== " << endl;
112
113	<	// Input fileName
114	<	TString infile1 = "skimmed361p3_20100623/Data_Incl_NewRelIso_wErrors_TandL_3";
115	<	TString infile2 = "skimmed361p3_20100623/MC_ppMuX_Incl_NewRelIso_wErrors_TandL_3";
116	<	TString infile3 = "skimmed361p3_20100623/MC_Wjets_Incl_NewRelIso_wErrors_TandL_3";
117	<	//   TString infile1 = "skimmed361p3_20100623/Data_Incl_OldRelIso_wErrors_TandL_3";
118	<	//   TString infile2 = "skimmed361p3_20100623/MC_ppMuX_Incl_OldRelIso_wErrors_TandL_3";
119	<	//   TString infile3 = "skimmed361p3_20100623/MC_Wjets_Incl_OldRelIso_wErrors_TandL_3";
113	>	// ////////////////////////
114	>	// Input/output fileName
115	>	// ////////////////////////
116	>	TString infile1,infile2,infile3,infile4,outFile;
117	>	if (isoname == "New") {//skimmed361p3_20100624: Loose ; skimmed361p3_20100623: tight
118	>	infile1 = "skimmed361p3_20100720/Data_250_Incl_NewRelIso_wErrors_TandL_4";
119	>	infile2 = "skimmed361p3_MC/MC_Mu15_Incl_NewRelIso_wErrors_TandL_4"; //InclusiveMu15
120	>	//infile2 = "skimmed361p3_MC/MC_ppMuX_Incl_NewRelIso_wErrors_TandL_4";
121	>	infile3 = "skimmed361p3_MC/MC_Wjets_Incl_NewRelIso_wErrors_TandL_4";
122	>	infile4 = "skimmed361p3_MC/MC_Zjets_Incl_NewRelIso_wErrors_TandL_4";
123	>
124	>	//outFile = "Plots/20100720/QCD_data_LandauFit_4_Incl_ppMuX_Scale0_250_125_"; //ppMuX
125	>	outFile = "Plots/20100720/QCD_data_LandauFit_4_Incl_Scale0_250_125_"; //InclusiveMu15
126	>	}
127	>	if (isoname == "Old") {
128	>	infile1 = "skimmed361p3_20100720/Data_250_Incl1_OldRelIso_wErrors_TandL_4";
129	>	infile2 = "skimmed361p3_MC/MC_Mu15_Incl1_OldRelIso_wErrors_TandL_4"; //InclusiveMu15
130	>	//infile2 = "skimmed361p3_MC/MC_ppMuX_Incl_OldRelIso_wErrors_TandL_4";
131	>	infile3 = "skimmed361p3_MC/MC_Wjets_Incl1_OldRelIso_wErrors_TandL_4";
132	>	infile4 = "skimmed361p3_MC/MC_Zjets_Incl1_OldRelIso_wErrors_TandL_4";
133
134	+	outFile = "Plots/20100720/QCD_data_GaussFit_4_Incl1_Scale0_250_";
135	+	}
136		// Output file
137	<	TString outFile = "Plots/20100623/QCD_data_LandauFit_3_";
113	<	//   TString outFile = "Plots/20100623/QCD_data_GaussFit_3_";
137	>
138		ofstream outTxtFile;
139		outTxtFile.open(TString(outFile+"Results.txt"));
140	+
141	+	double binwds[5] = {0.05,0.02,0.01,0.005,0.002};
142	+	double binwidth[4];
143	+	if (isoname == "New") {
144	+	for (int i=0; i<4;++i)
145	+	binwidth[i] = binwds[i];
146	+	}
147	+	else {
148	+	for (int i=0; i<4;++i)
149	+	binwidth[i] = binwds[i+1];
150	+	}
151
152	<	double binwidth[5] = {0.05,0.02,0.01,0.005,0.002};
153	<	int res_nbkg[5] = {0,0,0,0,0};
154	<	int res_nsig[5] = {0,0,0,0,0};
152	>	// ////////////////////////
153	>	// Looping of files of different binning
154	>	// ////////////////////////
155	>	int res_nbkg[4] = {0,0,0,0};
156	>	int res_nsig[4] = {0,0,0,0};
157		int nfiles = sizeof(binwidth)/sizeof(double);
158	+
159		for(int ibin = 0; ibin < nfiles ; ++ibin) {
160	+	// ////////////////////////
161	+	// define fit region:
162	+	// ////////////////////////
163	+	double ax0 = 0.0;
164	+	double axf = 0.0;
165	+	double bx0 = 0.0;
166	+	double bxf = 0.0;
167	+	double loosecut = 0.;
168	+
169	+	cout<<">>>>> Use "<<isoname<<" RelIso"<<endl;
170	+
171	+	if ( isoname.Contains("Old") ) {
172	+	// old
173	+	loosecut = 0.9; // 10/11
174	+	ax0 = 20./21.; // 0.95
175	+	axf = 1.0 + binwidth[ibin];
176	+	bx0 = 0.35; // Must be smaller than x at peak //0.35
177	+	bxf = 0.9; // Should not overlap signal region //0.9
178	+	}
179	+	else if ( isoname.Contains("New") ) {
180	+	// new
181	+	ax0 = 0.0;
182	+	axf = 0.05; // 1./19. = old 0.95
183	+	loosecut = 0.125; // 1./9. = old 0.9
184	+	bx0 = 0.1;//0.2
185	+	bxf = 1.;//1.2
186	+	}
187	+
188	+	//if (isoname == "Old" && ibin == 0) continue;
189		TString plotFile = outFile;
190	<	gStyle->SetOptFit(100);
190	>	if (drawStats) gStyle->SetOptFit(100);
191		cout << " =========================== " << endl;
192		cout << " = Bin width = " << binwidth[ibin] << endl;
193		cout << " =========================== " << endl;
#	Line 139 | Line 206 | void fitRelIso_Data()
206		// Data
207		m_fs["Data"] = pair<TString,double>(TString(infile1+"_"+suffix1.str()+".root"),1.);
208		// MC
209	<	m_fs["QCD"] = pair<TString,double>(TString(infile2+"_"+suffix1.str()+".root"),0.154166);
210	<	m_fs["WJets"] = pair<TString,double>(TString(infile3+"_"+suffix1.str()+".root"),0.0000369);
209	>	m_fs["QCD"] = pair<TString,double>(TString(infile2+"_"+suffix1.str()+".root"),0.0062348); //InclusiveMu15
210	>	//m_fs["QCD"] = pair<TString,double>(TString(infile2+"_"+suffix1.str()+".root"),2.9647348);
211	>	m_fs["WJets"] = pair<TString,double>(TString(infile3+"_"+suffix1.str()+".root"),0.0007934);
212	>	m_fs["ZJets"] = pair<TString,double>(TString(infile4+"_"+suffix1.str()+".root"),0.0007024);
213
214		// User fit function
215	<	TF1 *fuser = new TF1("fuser","[0]*TMath::Landau(x,[1],[2],0)",0.,2.);
216	<	fuser->SetParameters(1.,0.5,1.);
217	<	fuser->SetParLimits(2,0.125,1.);
218	<	fuser->SetParLimits(3,0.,0.5);
215	>	//     TF1 *fuser = new TF1("fuser","1./(1.+[0]*TMath::Landau(x,[1],[2],0))",ax0,axf);
216	>	//     TF1 *fu0 = new TF1("fu0","landau",0.,ax0);
217	>	//     TF1 *fu1 = new TF1("fu1","gaus",0.,ax0);
218	>	//     TF1 *fu2 = new TF1("fu2","landau",ax0,axf);
219	>	//     TF1 *fu3 = new TF1("fu3","landau(3)+gaus(2)",0.,ax0);
220	>	//     TF1 *fuser = new TF1("fuser","[0]*TMath::Landau(x,[1],[2],1)",ax0,axf);
221	>	TF1 *fuser = new TF1("fuser","landau",ax0,axf);
222	>	fuser->SetParameters(10.,0.3,1000.);//const ==> Remember to change according to L_{int}
223	>	fuser->SetParLimits(2,0.15,1.);//MPV
224	>	fuser->SetParLimits(3,0.,1.);//sigma
225	>	//fuser->SetParLimits(3,0.,0.5);
226
227		// File type
228		TString fileName = m_fs["Data"].first;
229		//cout << fileName << endl;
230
231		Int_t jbin;
232	<	Double_t na[3],na1[3];
232	>	Double_t na[4],na1[4];
233		int j = 0;
234		for (std::map<TString,pair<TString,double> >::iterator itFile = m_fs.begin();
235		itFile != m_fs.end(); ++itFile, ++j) {
#	Line 166 | Line 242 | void fitRelIso_Data()
242		tree->SetBranchAddress("na1",&na1[j]);
243		if ((*itFile).first == "Data") {
244		tree->SetBranchAddress("jbin",&jbin);
245	<	sf_mc = hs[(*itFile).first]->Integral();
245	>	if (isoname == "New") sf_mc = hs[(*itFile).first]->Integral(1,round(0.8/binwidth[ibin],0));
246	>	if (isoname == "Old") sf_mc = hs[(*itFile).first]->Integral(round(0.45/binwidth[ibin],0),hs[(*itFile).first]->GetMaximumBin());
247		//      sf_mc = hs[(*itFile).first]->GetEntries();
248	<	cout << ">>>>> Total events of data        = " << sf_mc << " <<<<<<<<<<<<<<<<" << endl;
248	>	if (scaleType == 0) cout << ">>>>> Total events of data        = " << sf_mc << " <<<<<<<<<<<<<<<<" << endl;
249		}
250		else
251		hs[(*itFile).first]->Scale(((*itFile).second).second);
#	Line 178 | Line 255 | void fitRelIso_Data()
255		// Calculate scale factor for MC; normalized to total number of events from data
256		hs["MCBKG"] = (TH1D*) hs["QCD"]->Clone();
257		hs["MCBKG"]->Add((TH1D*) hs["WJets"]->Clone());
258	<	double nTotMC = hs["MCBKG"]->Integral();
259	<	cout << ">>>>> Total weighted MC BKG       = " << nTotMC  << " <<<<<<<<<<<<<<<<" << endl;
258	>	hs["MCBKG"]->Add((TH1D*) hs["ZJets"]->Clone());
259	>	double nTotMC = 0;
260	>	if (isoname == "New") nTotMC = hs["MCBKG"]->Integral(1,round(0.8/binwidth[ibin],0));
261	>	if (isoname == "Old") nTotMC = hs["MCBKG"]->Integral(round(0.45/binwidth[ibin],0),hs["MCBKG"]->GetMaximumBin());
262	>	if (scaleType == 0) cout << ">>>>> Total weighted MC BKG       = " << nTotMC  << " <<<<<<<<<<<<<<<<" << endl;
263		//double nTotMC = hs["QCD"]->GetEntries()*m_fs["QCD"].second + hs["WJets"]->GetEntries()*m_fs["WJets"].second;
264		//cout << ">>>>> Total weighted MC BKG       = " << nTotMC << " <<<<<<<<<<<<<<<<" << endl;
265		sf_mc /= nTotMC;
266	<	cout << ">>>>> Scale factor for MC (sf_mc) = " << sf_mc << " <<<<<<<<<<<<<<<<" << endl;
187	<
188	<	hs["QCD"]->Scale(sf_mc);
189	<	hs["WJets"]->Scale(sf_mc);
190	<
191	<	// define fit region:
192	<	double ax0 = 0.0;
193	<	double axf = 0.0;
194	<	double bx0 = 0.0;
195	<	double bxf = 0.0;
196	<	double loosecut = 0.;
197	<
198	<	cout<<">>>>> Use "<<isoname<<" RelIso"<<endl;
199	<
200	<	if ( isoname.Contains("Old") ) {
201	<	// old
202	<	loosecut = 10./11.; // 0.9
203	<	ax0 = 20./21.; // 0.95
204	<	axf = 1.0 + binwidth[ibin];
205	<	bx0 = 0.3; // Must be smaller than x at peak
206	<	bxf = 0.9; // Should not overlap signal region
207	<	}
208	<	else if ( isoname.Contains("New") ) {
209	<	// new
210	<	ax0 = 0.0;
211	<	axf = 0.05; // 1./19. = old 0.95
212	<	loosecut = 0.125; // 1./9. = old 0.9
213	<	bx0 = 0.2;
214	<	bxf = 1.2;
215	<	}
266	>	if (scaleType == 0) cout << ">>>>> Scale factor for MC (sf_mc) = " << sf_mc << " <<<<<<<<<<<<<<<<" << endl;
267
268		int niter = 1;
269		double pass[3] = {0.,0.,0.};
#	Line 220 | Line 271 | void fitRelIso_Data()
271		double ns_tight[3] = {0.,0.,0.};
272		double ns_loose[3] = {0.,0.,0.};
273		double ns2 = 0.;
274	+	double chi2[2] = {0.,0.};
275	+	int ndof[2] = {0,0};
276
277		// Fit QCD in background region
225	–	cout<<">>>>> Fitting with "<<s0<<" function!"<<endl;
278		hs["data_all_clone"] = (TH1D*) hs["Data"]->Clone();
279	<	hs["data_all_clone"]->SetLineWidth(2);
279	>	//hs["data_all_clone"]->SetLineWidth(2);
280		hs["data_all_clone"]->SetMinimum(0);
281		hs["mc_qcd_clone"] = (TH1D*) hs["QCD"]->Clone();
282		hs["mc_wj_clone"] = (TH1D*) hs["WJets"]->Clone();
283	+	hs["mc_zj_clone"] = (TH1D*) hs["ZJets"]->Clone();
284
285		// Get num of bin at peak and bin width
286		Int_t nbxMax = hs["data_all_clone"]->GetXaxis()->FindBin(2.0);
287	<	hs["data_all_clone"]->GetXaxis()->SetRange(2,hs["data_all_clone"]->GetXaxis()->FindBin(bxf));
287	>	hs["data_all_clone"]->GetXaxis()->SetRange(round(0.05/binwidth[ibin],0),hs["data_all_clone"]->GetXaxis()->FindBin(bxf));
288		Int_t nbMax = hs["data_all_clone"]->GetMaximumBin();
289		cout << ">>>>> Bin number at peak = " << nbMax << endl;
290		hs["data_all_clone"]->GetXaxis()->SetRange(1,nbxMax);
#	Line 239 | Line 292 | void fitRelIso_Data()
292		double bwidth = axis0->GetBinWidth(nbMax);
293		cout << ">>>>> Histo bin width = " << bwidth << endl;
294
295	<	TString title = "RelIso";
296	<	if ( jbin == -1 )
297	<	title += " (inclusive)";
298	<	else if ( jbin == 4 )
299	<	title += " (#geq 4 jets)";
300	<	else if ( jbin == 1 ) {
301	<	title += " (";
302	<	title += jbin;
303	<	title += " jet)";
304	<	}
305	<	else if ( jbin < 4 ) {
306	<	title += " (";
307	<	title += jbin;
308	<	title += " jets)";
309	<	hs["data_all_clone"]->SetTitle(title);
310	<	}
311	<	else if ( jbin == 5 ) {
312	<	hs["data_all_clone"]->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
313	<	}
314	<	else {
315	<	cout<<"Don't mess around!!! jbin should be less than 6!"<<endl;
316	<	return;
317	<	}
295	>	TString title;
296	>	if (isoname == "New") title = "Relative Isolation";
297	>	if (isoname == "Old") title = "Relative Isolation (old)";
298	>	//     if ( jbin == -1 )
299	>	//       title += " (N_{jets} #geq 0)";
300	>	//     else if ( jbin == 5 ) {
301	>	//       title += "(#geq 1 jets (up to 3 jets))";
302	>	//     }
303	>	//     else if ( jbin > 10) {
304	>	//       ostringstream tmp_;
305	>	//       tmp_ << (jbin-10);
306	>	//       title += " (N_{jets} = ";
307	>	//       title += tmp_.str();
308	>	//       title += ")";
309	>	//     }
310	>	//     else if ( jbin <= 4 ) {
311	>	//       title += " (N_{jets} #geq ";
312	>	//       title += jbin;
313	>	//       title += ")";
314	>	//     }
315	>	//     else {
316	>	//       cout<<"jbin = " << jbin << " not supported!"<<endl;
317	>	//       return;
318	>	//     }
319
320	<	cvs[TString("cv"+suffix2.str())] = new TCanvas(TString("cv"+suffix2.str()),TString("cv"+suffix2.str()),600,600);
320	>	cvs[TString("cv"+suffix2.str())] = new TCanvas(TString("cv"+suffix2.str()),TString("cv"+suffix2.str()),10,10,700,500);
321		cvs[TString("cv"+suffix2.str())]->cd();
322
323		// Data
324		hs["data_all_clone"]->GetXaxis()->SetTitle(title);
325		hs["data_all_clone"]->GetXaxis()->SetLabelFont(42);
326	<	hs["data_all_clone"]->GetXaxis()->SetLabelSize(0.03);
326	>	hs["data_all_clone"]->GetXaxis()->SetLabelSize(0.04);
327		hs["data_all_clone"]->GetXaxis()->SetTitleFont(42);
328	<	hs["data_all_clone"]->GetXaxis()->SetTitleSize(0.035);
328	>	hs["data_all_clone"]->GetXaxis()->SetTitleSize(0.05);
329		hs["data_all_clone"]->GetXaxis()->SetTitleOffset(1.2);
330
331	<	hs["data_all_clone"]->GetYaxis()->SetTitle("Number of Events");
331	>	hs["data_all_clone"]->GetYaxis()->SetTitle(TString("Events / "+suffix1.str()));
332	>	//hs["data_all_clone"]->GetYaxis()->SetTitle(TString("Events"));
333		hs["data_all_clone"]->GetYaxis()->SetLabelFont(42);
334	<	hs["data_all_clone"]->GetYaxis()->SetLabelSize(0.03);
334	>	hs["data_all_clone"]->GetYaxis()->SetLabelSize(0.04);
335		hs["data_all_clone"]->GetYaxis()->SetTitleFont(42);
336	<	hs["data_all_clone"]->GetYaxis()->SetTitleSize(0.035);
336	>	hs["data_all_clone"]->GetYaxis()->SetTitleSize(0.05);
337		hs["data_all_clone"]->GetYaxis()->SetTitleOffset(1.4);
338	<	hs["data_all_clone"]->GetXaxis()->SetRangeUser(0.,plot_xMax);
338	>	hs["data_all_clone"]->GetXaxis()->SetRangeUser(plot_xMin,plot_xMax);
339		hs["data_all_clone"]->SetMarkerStyle(20);
340		hs["data_all_clone"]->SetMarkerSize(0.8);
341		hs["data_tmp_clone"] = (TH1D*) hs["data_all_clone"]->Clone(); // for zoomed in plot
342
288	–	if (DrawOverFlow) { // store histo with overflow bin drawn
289	–	int nbins_ = hs["data_all_clone"]->GetXaxis()->FindBin(plot_xMax);
290	–	for (std::map<TString,TH1*>::iterator itHS = hs.begin();
291	–	itHS != hs.end(); ++itHS) {
292	–	hso[(*itHS).first] = (TH1D*) ((*itHS).second)->Clone();
293	–	hso[(*itHS).first]->SetBins(nbins_,0.,plot_xMax+bwidth);
294	–	double v_ovf = 0.;
295	–	for (int i=0; i < ((*itHS).second)->GetNbinsX()+1 ; i++) {
296	–	if (i < nbins_) {
297	–	hso[(*itHS).first]->SetBinContent(i+1,hs[(*itHS).first]->GetBinContent(i+1));
298	–	hso[(*itHS).first]->SetBinError(i+1,hs[(*itHS).first]->GetBinError(i+1));
299	–	}
300	–	else {
301	–	v_ovf += hs[(*itHS).first]->GetBinContent(i+1);
302	–	}
303	–	}
304	–	hso[(*itHS).first]->SetBinContent(nbins_,v_ovf);
305	–	hso[(*itHS).first]->SetBinError(nbins_,TMath::Sqrt(v_ovf));
306	–	}
307	–	}
308	–
309	–	hst[TString("Data"+suffix2.str())] = new THStack(TString("Data"+suffix2.str()),TString("Data"+suffix2.str()+"stacked"));
310	–	if (DrawOverFlow) {
311	–	hso["mc_qcd_clone"]->SetStats(0);
312	–	hso["mc_qcd_clone"]->SetLineWidth(2);
313	–	hso["mc_qcd_clone"]->SetLineColor(40);
314	–	hso["mc_qcd_clone"]->SetFillStyle(3018);
315	–	hso["mc_qcd_clone"]->SetFillColor(38);
316	–
317	–	hso["mc_wj_clone"]->SetStats(0);
318	–	hso["mc_wj_clone"]->SetLineWidth(2);
319	–	hso["mc_wj_clone"]->SetLineColor(8);
320	–	//hso["mc_wj_clone"]->SetFillStyle(3001);
321	–	hso["mc_wj_clone"]->SetFillColor(3);
322	–	hst[TString("Data"+suffix2.str())]->Add(hso["mc_qcd_clone"]);
323	–	hst[TString("Data"+suffix2.str())]->Add(hso["mc_wj_clone"]);
324	–	}
325	–	else {
326	–	hs["mc_qcd_clone"]->SetStats(0);
327	–	hs["mc_qcd_clone"]->SetLineWidth(2);
328	–	hs["mc_qcd_clone"]->SetLineColor(40);
329	–	hs["mc_qcd_clone"]->SetFillStyle(3018);
330	–	hs["mc_qcd_clone"]->SetFillColor(38);
331	–
332	–	hs["mc_wj_clone"]->SetStats(0);
333	–	hs["mc_wj_clone"]->SetLineWidth(2);
334	–	hs["mc_wj_clone"]->SetLineColor(8);
335	–	//hs["mc_wj_clone"]->SetFillStyle(3001);
336	–	hs["mc_wj_clone"]->SetFillColor(3);
337	–
338	–	hst[TString("Data"+suffix2.str())]->Add(hs["mc_qcd_clone"]);
339	–	hst[TString("Data"+suffix2.str())]->Add(hs["mc_wj_clone"]);
340	–	}
341	–
343		// Find good y range first
344		if (isoname == "New") {
345		if ((hs["data_all_clone"]->GetBinContent(1) + hs["data_all_clone"]->GetBinError(1))*1.05 > yMax)
#	Line 363 | Line 364 | void fitRelIso_Data()
364		cout << ">>>>> Start fitting control region: " << bx0 << " to " << bxf << endl;
365
366		// Very first fit
367	+	cout << ">>>>> Fitting with " << s0 << " function!" << endl;
368	+	fuser->SetParameters(10.,0.3,100.);//const
369	+	fuser->SetParLimits(2,0.15,1.);//MPV
370	+	fuser->SetParLimits(3,0.,1.);//sigma
371		hs["data_all_clone"]->Fit(s0,fitOpt,"",bx0,bxf);
372
373		// Very first fit function
374		TF1 *f0 = (TF1*) hs["data_all_clone"]->GetFunction(s0)->Clone();
375		Int_t npar = f0->GetNpar();
376	<	Double_t chi2 = f0->GetChisquare();
377	<	Int_t ndof = f0->GetNDF();
378	<	pass[0] = chi2/ndof; // Very first fit norm chi2. Fixed range.
376	>	chi2[0] = f0->GetChisquare();
377	>	ndof[0] = f0->GetNDF();
378	>	pass[0] = chi2[0]/ndof[0]; // Very first fit norm chi2. Fixed range.
379		cout << ">>> Norm chi2 = " << pass[0] << endl;
380		if (!dynamic)
381		pass[1] = pass[0];
#	Line 395 | Line 400 | void fitRelIso_Data()
400		if(dynamic) {
401		bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2);
402		bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2);
403	+	cout << "[Before] bx0 = " << bx0 << endl;
404	+	cout << "[Before] bxf = " << bxf << endl;
405		}
406		if (isoname == "New") {
407		bx0 = ( bx0 <= loosecut ? loosecut : bx0 );
#	Line 407 | Line 414 | void fitRelIso_Data()
414		bxf = ( bxf > loosecut ? loosecut : bxf );
415		bx0 = ( bx0 < 0. ? 0. : bx0 );
416		}
417	+	//     if (isoname == "New") {
418	+	//       bx0 = ( bx0 <= loosecut ? loosecut :
419	+	//            ( bx0 < hs["data_all_clone"]->GetXaxis()->GetXmax() ? bx0 :
420	+	//              ( hs["data_all_clone"]->GetXaxis()->GetXmax() + loosecut )/2. ) );
421	+	//       bxf = ( bxf >= hs["data_all_clone"]->GetXaxis()->GetXmax() ?
422	+	//            hs["data_all_clone"]->GetXaxis()->GetXmax() :
423	+	//            ( bxf > bx0 ? bxf :
424	+	//              ( bx0 + hs["data_all_clone"]->GetXaxis()->GetXmax() )/2. ) );
425	+	//     }
426	+	//     if (isoname == "Old") {
427	+	//       bxf = ( bxf > loosecut ? loosecut :
428	+	//            ( bxf > hs["data_all_clone"]->GetXaxis()->GetXmin() ? bxf :
429	+	//              ( hs["data_all_clone"]->GetXaxis()->GetXmin() + loosecut)/2. ) );
430	+	//       bx0 = ( bx0 < hs["data_all_clone"]->GetXaxis()->GetXmin() ? hs["data_all_clone"]->GetXaxis()->GetXmin() :
431	+	//            ( bx0 > bxf ? ( bxf + hs["data_all_clone"]->GetXaxis()->GetXmin() )/2. : bx0) );
432	+	//     }
433		bound[0] = bx0;
434		bound[1] = bxf;
435		if (dynamic) {
#	Line 416 | Line 439 | void fitRelIso_Data()
439
440		double delta = pass[1]-pass[2];
441
442	<	while (dynamic && niter <= 20 && abs(delta) > 0.000001 && ndof >= 4) {
442	>	while (dynamic && niter <= 20 && abs(delta) > 0.000001 && ndof[0] >= 4) {
443		if (niter > 1)
444		pass[1] = pass[2];
445		if (delta != 0){
#	Line 431 | Line 454 | void fitRelIso_Data()
454		cout << "\n>>> Start fitting new range (bx0, bxf) = (" << bx0
455		<< ", " << bxf << ")" << endl;
456
457	+	fuser->SetParameters(10.,0.3,100.);//const
458	+	fuser->SetParLimits(2,0.15,1.);//MPV
459	+	fuser->SetParLimits(3,0.,1.);//sigma
460	+
461		hs["data_all_clone"]->Fit(s0,fitOpt,"",bx0,bxf);
462		TF1 *fi = (TF1*) hs["data_all_clone"]->GetFunction(s0)->Clone();
463		cout << " >> fi->GetChisquare() = " << fi->GetChisquare() << "; fi->GetNDF() = " << fi->GetNDF() << endl;
464	<	if ( fi->GetNDF() > 0 ) pass[2] = fi->GetChisquare()/fi->GetNDF();
465	<	else pass[2] = 999.;
464	>	if ( fi->GetNDF() > 0 )
465	>	pass[2] = fi->GetChisquare()/fi->GetNDF();
466	>	else {
467	>	chi2[1] = 999.;
468	>	ndof[1] = 1.;
469	>	pass[2] = 999.;
470	>	}
471		cout << " >> Current step norm. chi2  = " << pass[2] << endl;
472		delta = pass[1]-pass[2];
473		cout << " >> Delta = " << delta << endl;
474
475		if (delta > 0 && bx0 >= loosecut && bxf <= hs["data_all_clone"]->GetXaxis()->GetXmax() && bx0 < bxf) {
476	+	chi2[1] = fi->GetChisquare();
477	+	ndof[1] = fi->GetNDF();
478		bound[0] = bx0;
479		bound[1] = bxf;
480		printf(" >> Function has %i parameters. Chisquare = %g\n",
#	Line 458 | Line 492 | void fitRelIso_Data()
492		}
493		bx0 = fi->GetParameter(1) - wxmin*fi->GetParameter(2);
494		bxf = fi->GetParameter(1) + wxmax*fi->GetParameter(2);
495	+	cout << "[Before] bx0 = " << bx0 << endl;
496	+	cout << "[Before] bxf = " << bxf << endl;
497	+
498		if (isoname == "New") {
499		bx0 = ( bx0 < loosecut ? loosecut : bx0 );
500		bxf = ( bxf > 2. ? 2. : ( bxf > bx0 ? bxf : 1.1*bx0 ) );
#	Line 466 | Line 503 | void fitRelIso_Data()
503		bxf = ( bxf > loosecut ? loosecut : bxf );
504		bxf = ( bx0 < 0. ? 0. : bx0 );
505		}
506	+	//        if (isoname == "New") {
507	+	//          bx0 = ( bx0 <= loosecut ? loosecut :
508	+	//                  ( bx0 < hs["data_all_clone"]->GetXaxis()->GetXmax() ? bx0 :
509	+	//                    ( hs["data_all_clone"]->GetXaxis()->GetXmax() + loosecut )/2. ) );
510	+	//          bxf = ( bxf >= hs["data_all_clone"]->GetXaxis()->GetXmax() ?
511	+	//                  hs["data_all_clone"]->GetXaxis()->GetXmax() :
512	+	//                  ( bxf > bx0 ? bxf :
513	+	//                    ( bx0 + hs["data_all_clone"]->GetXaxis()->GetXmax() )/2. ) );
514	+	//        }
515	+	//        if (isoname == "Old") {
516	+	//          bxf = ( bxf > loosecut ? loosecut :
517	+	//                  ( bxf > hs["data_all_clone"]->GetXaxis()->GetXmin() ? bxf :
518	+	//                    ( hs["data_all_clone"]->GetXaxis()->GetXmin() + loosecut)/2. ) );
519	+	//          bx0 = ( bx0 < hs["data_all_clone"]->GetXaxis()->GetXmin() ? hs["data_all_clone"]->GetXaxis()->GetXmin() :
520	+	//                  ( bx0 > bxf ? ( bxf + hs["data_all_clone"]->GetXaxis()->GetXmin() )/2. :bx0 ) );
521	+	//        }
522		cout << ">>> Range for next iteration (bx0, bxf) = (" << bx0;
523		cout << ", " << bxf << ")" << endl;
524		}
#	Line 482 | Line 535 | void fitRelIso_Data()
535		pass[0]=pass[2]; // First dynamic fit norm. chi2
536		} // delta != 0
537		} // End while
538	+
539	+	// Output equivalent range
540	+	if (dynamic) {
541	+	if (bound[0]/binwidth[ibin] - floor(bound[0]/binwidth[ibin]) > 0.5)
542	+	bound[0] = round(bound[0]/binwidth[ibin],0) * binwidth[ibin];
543	+	else
544	+	bound[0] = (floor(bound[0]/binwidth[ibin])+0.5) * binwidth[ibin];
545	+	if (bound[1]/binwidth[ibin] - floor(bound[1]/binwidth[ibin]) >= 0.5)
546	+	bound[1] = round(bound[1]/binwidth[ibin],0) * binwidth[ibin];
547	+	else
548	+	bound[1] = round(bound[1]/binwidth[ibin],0) * binwidth[ibin];
549	+	}
550	+
551		if(debug){
552		cout << ">>> First dynamic fit norm. chi2   = " << pass[0] << endl;
553		cout << ">>> Min dynamic norm. chi2         = " << pass[1] << endl;
#	Line 500 | Line 566 | void fitRelIso_Data()
566		nfac1 -= (hs["data_all_clone"]->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bound[1])/axis->GetBinWidth(bmax);
567		cout << ">>> Final Nfac = " << nfac1 << endl;
568
503	–	// ////////////////////////
504	–	// Histos and extrapolation
505	–	// ////////////////////////
506	–
507	–	if (DrawOverFlow) {
508	–	hso["data_all_clone"]->GetYaxis()->SetRangeUser(0.,yMax);
509	–	hso["data_all_clone"]->Draw("esame");
510	–	}
511	–
512	–	hst[TString("Data"+suffix2.str())]->Draw("sames hist");
513	–
569		TF1 *f1 = (TF1*) f0->Clone();
570		f1->SetRange(bound[0],bound[1]);
571
#	Line 521 | Line 576 | void fitRelIso_Data()
576		}
577		f1->SetLineColor(2);
578
524	–	// Draw on top of everything
525	–	gStyle->SetErrorX(0.5);
526	–	if (DrawOverFlow) {
527	–	hso["data_all_clone"]->Draw("esame");
528	–	}
529	–	else {
530	–	hs["data_all_clone"]->Draw("esame"); // Draw on top of others
531	–	}
532	–	//f1->SetLineWidth(3.5);
533	–	f1->Draw("same");
534	–
535	–	// Connect fitted and extrapolated region
536	–	TF1 *fin = (TF1*) f1->Clone();
537	–	if (isoname == "New") {
538	–	fin->SetRange(axf,bound[0]);
539	–	fin->SetLineStyle(2);
540	–	fin->SetLineColor(9);
541	–	fin->SetLineWidth(3.5);
542	–	fin->Draw("same");
543	–
544	–	if (drawtail){
545	–	TF1 *fine = (TF1*) f1->Clone();
546	–	fine->SetRange(bound[1],2.0);
547	–	fine->SetLineStyle(2);
548	–	fine->SetLineColor(9);
549	–	fine->SetLineWidth(3.5);
550	–	fine->Draw("same");
551	–	}
552	–	}
553	–	if (isoname == "Old") {
554	–	fin->SetRange(bound[1],axf);
555	–	fin->SetLineStyle(2);
556	–	fin->SetLineColor(9);
557	–	fin->SetLineWidth(3.5);
558	–	fin->Draw("same");
559	–
560	–	if (drawtail){
561	–	TF1 *fine = (TF1*) f1->Clone();
562	–	fine->SetRange(0.,bound[0]);
563	–	fine->SetLineStyle(2);
564	–	fine->SetLineColor(9);
565	–	fine->SetLineWidth(3.5);
566	–	fine->Draw("same");
567	–	}
568	–	}
569	–
579		TF1 *f2 = (TF1*) f1->Clone();
580		f2->SetRange(ax0,axf);
581		f2->SetLineColor(4);
582		//f2->SetLineWidth(3.5);
574	–	f2->Draw("same");
583
584		// ////////////////////////
585		// Extract number of qcd events
#	Line 634 | Line 642 | void fitRelIso_Data()
642		* fabs(axf-ax0)/(2*bwidth);
643		ns_loose[1] = (TMath::Landau(ax0,par0[1],par0[2],kFALSE)*par0[0] + TMath::Landau(loosecut,par0[1],par0[2],kFALSE)*par0[0])
644		* fabs(loosecut-ax0)/(2*bwidth);
645	<	// Method 3: area of Landau function
645	>	// Method 3: area of Landau function (default)
646		ns_tight[2] = f2->Integral(ax0,axf)/bwidth;
647		ns_loose[2] = f2->Integral(ax0,loosecut)/bwidth;
648
649	<	// Stats box
650	<	TPaveStats * tps1;
651	<	if(!dynamic){
652	<	cvs[TString("cv"+suffix2.str())]->Update();
653	<	tps1 = (TPaveStats*) cvs[TString("cv"+suffix2.str())]->GetPrimitive("stats");
649	>	// ////////////////////////
650	>	// Histos and extrapolation
651	>	// ////////////////////////
652	>
653	>	// Draw stacked MC
654	>	double sf_mc1 = 0.;
655	>	if (isoname == "New") sf_mc1 = ns_tight[2]/hs["QCD"]->Integral(1,round(0.05/binwidth[ibin],0));
656	>	if (isoname == "Old") sf_mc1 = ns_tight[2]/hs["QCD"]->Integral(round(0.95/binwidth[ibin],0)+1,hs["QCD"]->GetMaximumBin());
657	>	if (scaleType == 0) {
658	>	hs["mc_qcd_clone"]->Scale(sf_mc);
659	>	hs["mc_wj_clone"]->Scale(sf_mc);
660	>	}
661	>	if (scaleType == 1) {
662	>	sf_mc = sf_mc1;
663	>	hs["mc_qcd_clone"]->Scale(sf_mc1);
664	>	hs["mc_wj_clone"]->Scale(sf_mc1);
665	>	cout << ">>>>> Scale factor for MC (sf_mc) = " << sf_mc1 << " <<<<<<<<<<<<<<<<" << endl;
666	>	}
667	>
668	>	if (DrawOverFlow) { // store histo with overflow bin drawn
669	>	int nbins_ = hs["data_all_clone"]->GetXaxis()->FindBin(plot_xMax);
670	>	for (std::map<TString,TH1*>::iterator itHS = hs.begin();
671	>	itHS != hs.end(); ++itHS) {
672	>	hso[(*itHS).first] = (TH1D*) ((*itHS).second)->Clone();
673	>	hso[(*itHS).first]->SetBins(nbins_,0.,plot_xMax+bwidth);
674	>	double v_ovf = 0.;
675	>	for (int i=0; i < ((*itHS).second)->GetNbinsX()+1 ; i++) {
676	>	if (i < nbins_) {
677	>	hso[(*itHS).first]->SetBinContent(i+1,hs[(*itHS).first]->GetBinContent(i+1));
678	>	hso[(*itHS).first]->SetBinError(i+1,hs[(*itHS).first]->GetBinError(i+1));
679	>	}
680	>	else {
681	>	v_ovf += hs[(*itHS).first]->GetBinContent(i+1);
682	>	}
683	>	}
684	>	hso[(*itHS).first]->SetBinContent(nbins_,v_ovf);
685	>	hso[(*itHS).first]->SetBinError(nbins_,TMath::Sqrt(v_ovf));
686	>	}
687	>	}
688	>
689	>	hst[TString("Data"+suffix2.str())] = new THStack(TString("Data"+suffix2.str()),TString("Data"+suffix2.str()+"stacked"));
690	>	if (DrawOverFlow) {
691	>	hso["mc_qcd_clone"]->SetStats(0);
692	>	//hso["mc_qcd_clone"]->SetLineWidth(2);
693	>	//hso["mc_qcd_clone"]->SetLineColor(style.QCDColor);
694	>	hso["mc_qcd_clone"]->SetFillColor(style.QCDColor);
695	>	hso["mc_qcd_clone"]->SetFillStyle(style.QCDFill);
696	>
697	>	hso["mc_zj_clone"]->SetStats(0);
698	>	hso["mc_zj_clone"]->SetFillStyle(style.DYZJetsFill);
699	>	hso["mc_zj_clone"]->SetFillColor(style.DYZJetsColor);
700	>
701	>	hso["mc_wj_clone"]->SetStats(0);
702	>	//hso["mc_wj_clone"]->SetLineWidth(2);
703	>	//hso["mc_wj_clone"]->SetLineColor(style.WJetsColor);
704	>	hso["mc_wj_clone"]->SetFillStyle(style.WJetsFill);
705	>	hso["mc_wj_clone"]->SetFillColor(style.WJetsColor);
706	>	hst[TString("Data"+suffix2.str())]->Add(hso["mc_qcd_clone"]);
707	>	hst[TString("Data"+suffix2.str())]->Add(hso["mc_zj_clone"]);
708	>	hst[TString("Data"+suffix2.str())]->Add(hso["mc_wj_clone"]);
709		}
710		else {
711	<	hs["data_all_clone"]->Fit(s0,fitOpt,"",bound[0],bound[1]);
712	<	cvs[TString("cv"+suffix2.str())]->Update();
713	<	tps1 = (TPaveStats*) cvs[TString("cv"+suffix2.str())]->GetPrimitive("stats");
714	<	}
715	<
716	<	tps1->SetName(TString("Landau fit"+suffix2.str()));
717	<
718	<
719	<	TList *list = tps1->GetListOfLines();
720	<	//   TText *tconst = tps1->GetLineWith("Constant");
721	<	//   list->Remove(tconst);
722	<	TLatex *myt = new TLatex(0,0,"Landau Fit");
723	<	//TLatex *myt = new TLatex(0,0,"#font[22]{Landau Fit}");
724	<	list->AddFirst(myt);
725	<	hs["data_all_clone"]->SetStats(0);
726	<	gStyle->SetOptFit(0);
727	<	cvs[TString("cv"+suffix2.str())]->Modified();
728	<
729	<	tps1->SetTextFont(42);
730	<	tps1->SetTextSize(0.03);
731	<	tps1->SetFillColor(0);
732	<	tps1->SetFillStyle(0);
733	<	tps1->SetBorderSize(0);
734	<	tps1->SetX1NDC(0.63);
735	<	tps1->SetX2NDC(0.88);
736	<	tps1->SetY1NDC(0.74);
737	<	tps1->SetY2NDC(0.87);
738	<	cvs[TString("cv"+suffix2.str())]->cd();
739	<	gPad->Update();
711	>	hs["mc_qcd_clone"]->SetStats(0);
712	>	//hs["mc_qcd_clone"]->SetLineWidth(2);
713	>	//hs["mc_qcd_clone"]->SetLineColor(style.QCDColor);
714	>	hs["mc_qcd_clone"]->SetFillStyle(style.QCDFill);
715	>	hs["mc_qcd_clone"]->SetFillColor(style.QCDColor);
716	>
717	>	hs["mc_zj_clone"]->SetStats(0);
718	>	hs["mc_zj_clone"]->SetFillStyle(style.DYZJetsFill);
719	>	hs["mc_zj_clone"]->SetFillColor(style.DYZJetsColor);
720	>
721	>	hs["mc_wj_clone"]->SetStats(0);
722	>	//hs["mc_wj_clone"]->SetLineWidth(2);
723	>	//hs["mc_wj_clone"]->SetLineColor(style.WJetsColor);
724	>	hs["mc_wj_clone"]->SetFillStyle(style.WJetsFill);
725	>	hs["mc_wj_clone"]->SetFillColor(style.WJetsColor);
726	>
727	>	hst[TString("Data"+suffix2.str())]->Add(hs["mc_qcd_clone"]);
728	>	hst[TString("Data"+suffix2.str())]->Add(hs["mc_zj_clone"]);
729	>	hst[TString("Data"+suffix2.str())]->Add(hs["mc_wj_clone"]);
730	>	}
731	>
732	>	// Draw data
733	>	if (DrawOverFlow) {
734	>	hso["data_all_clone"]->GetYaxis()->SetRangeUser(0.,yMax);
735	>	hso["data_all_clone"]->Draw("esame");
736	>	}
737	>
738	>	hst[TString("Data"+suffix2.str())]->Draw("sames hist");
739	>
740	>	// Draw data again
741	>	gStyle->SetErrorX(0.5);
742	>
743	>	if (DrawOverFlow) {
744	>	hso["data_all_clone"]->Draw("esame");
745	>	}
746	>	else {
747	>	hs["data_all_clone"]->DrawClone("esame"); // Draw on top of others
748	>	}
749	>
750	>	// Draw fit function
751	>
752	>	// fit function in control region
753	>	//f1->SetLineWidth(3.5);
754	>	f1->Draw("same");
755	>
756	>	// fit fuction in signal region
757	>	f2->Draw("same");
758	>
759	>	// Connect fit function between regions and tail
760	>	TF1 *fin = (TF1*) f1->Clone();
761	>	if (isoname == "New") {
762	>	fin->SetRange(axf,bound[0]);
763	>	fin->SetLineStyle(2);
764	>	fin->SetLineColor(9);
765	>	fin->SetLineWidth(3.5);
766	>	fin->Draw("same");
767	>
768	>	if (drawtail){
769	>	TF1 *fine = (TF1*) f1->Clone();
770	>	fine->SetRange(bound[1],2.0);
771	>	fine->SetLineStyle(2);
772	>	fine->SetLineColor(9);
773	>	fine->SetLineWidth(3.5);
774	>	fine->Draw("same");
775	>	}
776	>	}
777	>	if (isoname == "Old") {
778	>	fin->SetRange(bound[1],ax0);
779	>	fin->SetLineStyle(2);
780	>	fin->SetLineColor(9);
781	>	fin->SetLineWidth(3.5);
782	>	fin->Draw("same");
783	>
784	>	if (drawtail){
785	>	TF1 *fine = (TF1*) f1->Clone();
786	>	fine->SetRange(0.,bound[0]);
787	>	fine->SetLineStyle(2);
788	>	fine->SetLineColor(9);
789	>	fine->SetLineWidth(3.5);
790	>	fine->Draw("same");
791	>	}
792	>	}
793	>
794	>	TString mytxt;
795	>	if (isoname == "New") mytxt = "Landau Fit";
796	>	if (isoname == "Old") mytxt = "Gaussian Fit";
797	>
798	>	if (drawStats) {
799	>	// Draw Stats box
800	>	TPaveStats * tps1;
801	>	if(!dynamic){
802	>	cvs[TString("cv"+suffix2.str())]->Update();
803	>	tps1 = (TPaveStats*) cvs[TString("cv"+suffix2.str())]->GetPrimitive("stats");
804	>	}
805	>	else {
806	>	hs["data_all_clone"]->Fit(s0,fitOpt,"",bound[0],bound[1]);
807	>	cvs[TString("cv"+suffix2.str())]->Update();
808	>	tps1 = (TPaveStats*) cvs[TString("cv"+suffix2.str())]->GetPrimitive("stats");
809	>	}
810	>
811	>	tps1->SetName(TString(s0+" fit"+suffix2.str()));
812
813	<	//cmsPrel(13);
813	>	TList *list = tps1->GetListOfLines();
814	>	//   TText *tconst = tps1->GetLineWith("Constant");
815	>	//   list->Remove(tconst);
816	>	TLatex *myt = new TLatex(0,0,mytxt);
817	>	//TLatex *myt = new TLatex(0,0,"#font[22]{Landau Fit}");
818	>	list->AddFirst(myt);
819	>	hs["data_all_clone"]->SetStats(0);
820	>	gStyle->SetOptFit(0);
821	>	cvs[TString("cv"+suffix2.str())]->Modified();
822	>
823	>	tps1->SetTextFont(42);
824	>	tps1->SetTextSize(0.03);
825	>	tps1->SetFillColor(0);
826	>	tps1->SetFillStyle(0);
827	>	tps1->SetBorderSize(0);
828	>	tps1->SetX1NDC(0.63);
829	>	tps1->SetX2NDC(0.88);
830	>	tps1->SetY1NDC(0.74);
831	>	tps1->SetY2NDC(0.87);
832	>	cvs[TString("cv"+suffix2.str())]->cd();
833	>	gPad->Update();
834	>	}
835	>
836	>	TLatex *text1, *text2;
837	>	text1 = new TLatex(3.570061, 23.08044, "CMS Preliminary");
838	>	text1->SetNDC();
839	>	text1->SetTextAlign(13);
840	>	if (isoname == "New")
841	>	text1->SetX(0.24);
842	>	else
843	>	text1->SetX(0.54);
844	>	text1->SetY(0.92);
845	>	//text1->SetLineWidth(2);
846	>	text1->SetTextFont(42);
847	>	text1->SetTextSize(0.05);
848	>	text1->SetTextSizePixels(24);// dflt=28
849	>	text1->Draw();
850	>
851	>	text2 = new TLatex(13.570061, 23.08044, "0.25 pb^{-1} at #sqrt{s} = 7 TeV");
852	>	text2->SetNDC();
853	>	text2->SetTextAlign(13);
854	>	if (isoname == "New")
855	>	text2->SetX(0.24);
856	>	else
857	>	text2->SetX(0.54);
858	>	text2->SetY(0.85);
859	>	//text2->SetLineWidth(2);
860	>	text2->SetTextFont(42);
861	>	text2->SetTextSize(0.05);
862	>	text2->SetTextSizePixels(24);// dflt=28
863	>	text2->Draw();
864	>
865	>	//cmsPrel(35);
866		// Label histo
867	<	TLegend * leg1 = new TLegend(0.18,0.47,0.45,0.87);
868	<	leg1->SetHeader("#font[22]{CMS preliminary 2010   #sqrt{s} = 7 TeV }");
867	>	//TLegend * leg1 = new TLegend(0.58,0.45,0.85,0.8);
868	>	TLegend * leg1;
869	>	if (isoname == "New")
870	>	leg1 = new TLegend(0.66,0.52,0.95,0.92);
871	>	else
872	>	leg1 = new TLegend(0.28,0.57,0.55,0.92);
873	>	//leg1->SetHeader("#font[22]{CMS preliminary 2010   #sqrt{s} = 7 TeV }");
874		leg1->SetTextSize(0.03);
875		leg1->SetFillColor(0);
876		leg1->SetFillStyle(0);
877		if (DrawOverFlow) {
878	<	leg1->AddEntry(hso["data_all_clone"],"Data   #int#font[12]{L}dt = 13 nb^{-1}","p");
879	<	leg1->AddEntry(f1,"Fit of all events in control region","l");
880	<	leg1->AddEntry(f2,"Extrapolation to signal region","l");
881	<	leg1->AddEntry(hso["mc_qcd_clone"],"QCD","f");
882	<	leg1->AddEntry(hso["mc_wj_clone"],"W+Jets","f");
878	>	//leg1->AddEntry(hso["data_all_clone"],"Data   #int#font[12]{L}dt = 35 nb^{-1}","p");
879	>	leg1->AddEntry(hso["data_all_clone"],"Data","lp");
880	>	//       leg1->AddEntry(f1,TString(mytxt+" in control region"),"l");
881	>	//       leg1->AddEntry(f2,"Extrapolation to signal region","l");
882	>	leg1->AddEntry(f1,mytxt,"l");
883	>	leg1->AddEntry(f2,"Extrapolation","l");
884	>	leg1->AddEntry(hso["mc_wj_clone"],style.WJetsText,"f");
885	>	leg1->AddEntry(hso["mc_zj_clone"],style.DYZJetsText,"f");
886	>	leg1->AddEntry(hso["mc_qcd_clone"],style.QCDText,"f");
887		}
888		else {
889	<	leg1->AddEntry(hs["data_all_clone"],"Data   #int#font[12]{L}dt = 13 nb^{-1}","p");
890	<	leg1->AddEntry(f1,"Fit of all events in control region","l");
891	<	leg1->AddEntry(f2,"Extrapolation to signal region","l");
892	<	leg1->AddEntry(hs["mc_qcd_clone"],"QCD","f");
893	<	leg1->AddEntry(hs["mc_wj_clone"],"W+Jets","f");
889	>	//leg1->AddEntry(hs["data_all_clone"],"Data   #int#font[12]{L}dt = 35 nb^{-1}","p");
890	>	leg1->AddEntry(hs["data_all_clone"],"Data","lp");
891	>	leg1->AddEntry(f1,mytxt,"l");
892	>	leg1->AddEntry(f2,"Extrapolation","l");
893	>	leg1->AddEntry(hs["mc_wj_clone"],style.WJetsText,"f");
894	>	leg1->AddEntry(hs["mc_zj_clone"],style.DYZJetsText,"f");
895	>	leg1->AddEntry(hs["mc_qcd_clone"],style.QCDText,"f");
896		}
897		leg1->Draw();
898
899		plotFile = plotFile + suffix1.str() + ".pdf";
900	<	//cvs[TString("cv"+suffix2.str())]->Print(plotFile);
900	>	//     if (ibin == 1 || ibin == 2)
901	>	cvs[TString("cv"+suffix2.str())]->Print(plotFile);
902
903		// Print results
904		cout << "\n";
#	Line 736 | Line 935 | void fitRelIso_Data()
935		//     cout << ">>>>> Total expected events Loose Na  = " << na1[0] << endl;
936		//     cout << ">>>>> Num of Signal events            = " << na1[0] - round(ns_loose[2],0) << endl;
937		//   }
739	–
938		//   cout << "Landau(0.05) = " << TMath::Landau(0.05,par0[1],par0[2],kFALSE)*par0[0] << endl;
939		//   cout << "Landau(0.1) = " << TMath::Landau(0.1,par0[1],par0[2],kFALSE)*par0[0] << endl;
940		cout << " ===========================\n\n " << endl;
941
942		outTxtFile << "===========================" << endl;
943	<	outTxtFile << ">>>>> Bin width = " << bwidth << endl;
944	<	//outTxtFile << ">>>>> Calculated Ns (ratio)     = " << ns_tight[0] << endl;
945	<	//outTxtFile << ">>>>> Calculated Ns (trapezoid) = " << ns_tight[1] << endl;
946	<	outTxtFile << ">>>>> Total number of events    = " << na[0] << endl;
947	<	outTxtFile << ">>>>> Calculated bkg (integral) = " << round(ns_tight[2],0) << endl;
948	<	outTxtFile << ">>>>> Num of Signal events      = " << na[0] - round(ns_tight[2],0) << endl;
949	<	outTxtFile << ">>>>> Total MC QCD expected     = " << round(na[1]*((m_fs["QCD"]).second)*sf_mc,0) << endl;
950	<	outTxtFile << ">>>>> Total MC WJets expected   = " << round(na[2]*((m_fs["WJets"]).second)*sf_mc,0) << endl;
943	>	outTxtFile << ">>>>> Bin width                  = " << bwidth << endl;
944	>	//outTxtFile << ">>>>> Calculated Ns (ratio)      = " << ns_tight[0] << endl;
945	>	//outTxtFile << ">>>>> Calculated Ns (trapezoid)  = " << ns_tight[1] << endl;
946	>	if (dynamic) {
947	>	outTxtFile << ">>>>> Min dynamic norm. chi2     = " << pass[1] << endl;
948	>	outTxtFile << ">>>>> chi2/ndof                  = " << chi2[1] << "/" << ndof[1] << " = "
949	>	<< chi2[1]/ndof[1] << endl;
950	>	}
951	>	else {
952	>	outTxtFile << ">>>>> Min norm. chi2             = " << pass[1] << endl;
953	>	outTxtFile << ">>>>> chi2/ndof                  = " << chi2[0] << "/" << ndof[0] << " = " << chi2[0]/ndof[0] << endl;
954	>	}
955	>	outTxtFile << ">>>>> Best fitting region        = (" << bound[0];
956	>	outTxtFile << " - " << bound[1] << ")" << endl;
957	>	outTxtFile << ">>>>> Total number of events     = " << na[0] << endl;
958	>	outTxtFile << ">>>>> Calculated bkg (integral)  = " << round(ns_tight[2],0) << endl;
959	>	outTxtFile << ">>>>> Num of Signal events       = " << na[0] - round(ns_tight[2],0) << endl;
960	>	outTxtFile << ">>>>> Total MC QCD expected(N)   = " << round(na[1]*((m_fs["QCD"]).second)*sf_mc,1) << " +/- " <<
961	>	round(TMath::Sqrt(na[1])*((m_fs["QCD"]).second)*sf_mc,1) << endl;
962	>	outTxtFile << ">>>>> Total MC WJets expected(N) = " << round(na[2]*((m_fs["WJets"]).second)*sf_mc,1) << " +/- " <<
963	>	round(TMath::Sqrt(na[2])*((m_fs["WJets"]).second)*sf_mc,1) << endl;
964	>	outTxtFile << ">>>>> Total MC ZJets expected(N) = " << round(na[3]*((m_fs["ZJets"]).second)*sf_mc,1) << " +/- " <<
965	>	round(TMath::Sqrt(na[3])*((m_fs["ZJets"]).second)*sf_mc,1) << endl;
966	>	outTxtFile << ">>>>> Total MC QCD expected(R)   = " << round(na[1]*((m_fs["QCD"]).second),1) << " +/- " <<
967	>	round(TMath::Sqrt(na[1])*((m_fs["QCD"]).second),1)<< endl;
968	>	outTxtFile << ">>>>> Total MC WJets expected(R) = " << round(na[2]*((m_fs["WJets"]).second),1) << " +/- " <<
969	>	round(TMath::Sqrt(na[2])*((m_fs["WJets"]).second),1) << endl;
970	>	outTxtFile << ">>>>> Total MC ZJets expected(R) = " << round(na[3]*((m_fs["ZJets"]).second),1) << " +/- " <<
971	>	round(TMath::Sqrt(na[3])*((m_fs["ZJets"]).second),2) << endl;
972		outTxtFile << "===========================\n\n " << endl;
973
974
975	<	cvs[TString("cv1"+suffix2.str())] = new TCanvas(TString("cv1"+suffix2.str()),TString("cv1"+suffix2.str()),600,600);
975	>	// /////////////
976	>	// Zoom-in plot
977	>	// /////////////
978	>	cvs[TString("cv1"+suffix2.str())] = new TCanvas(TString("cv1"+suffix2.str()),TString("cv1"+suffix2.str()),10,10,700,500);
979		cvs[TString("cv1"+suffix2.str())]->cd();
980	<	if (isoname == "New")
981	<	hs["data_tmp_clone"]->GetXaxis()->SetRangeUser(0.,0.1);
982	<	if (isoname == "Old")
980	>
981	>	if (isoname == "New") {
982	>	hs["data_tmp_clone"]->GetXaxis()->SetRangeUser(0.,0.06);
983	>	hs["data_tmp_clone"]->GetYaxis()->SetRangeUser(0.3,10000.*binwidth[ibin]/binwidth[3]);
984	>	}
985	>	if (isoname == "Old") {
986		hs["data_tmp_clone"]->GetXaxis()->SetRangeUser(0.9,axf);
987	<	hs["data_tmp_clone"]->GetYaxis()->SetRangeUser(0.,yMax);
987	>	hs["data_tmp_clone"]->GetYaxis()->SetRangeUser(0.3,10000.*binwidth[ibin]/binwidth[3]);
988	>	}
989		hs["data_tmp_clone"]->Draw("e");
990		hst[TString("Data"+suffix2.str())]->Draw("sames hist");
991		hs["data_tmp_clone"]->Draw("esame");
992
993	<	TLegend * leg2 = new TLegend(0.38,0.6,0.65,0.85);
994	<	leg2->SetHeader("#font[22]{CMS preliminary 2010   #sqrt{s} = 7 TeV }");
993	>	text1->Draw();
994	>	text2->Draw();
995	>
996	>	TLegend * leg2;
997	>	if (isoname == "New")
998	>	leg2 = new TLegend(0.66,0.52,0.95,0.92);
999	>	else
1000	>	leg2 = new TLegend(0.28,0.57,0.55,0.92);
1001	>	//leg2->SetHeader("#font[22]{CMS preliminary 2010   #sqrt{s} = 7 TeV }");
1002		leg2->SetTextSize(0.03);
1003		leg2->SetFillColor(0);
1004		leg2->SetFillStyle(0);
1005	+	//leg2->AddEntry(hs["data_all_clone"],"Data   #int#font[12]{L}dt = 35 nb^{-1}","p");
1006		if (DrawOverFlow) {
1007	<	leg2->AddEntry(hso["data_all_clone"],"Data   #int#font[12]{L}dt = 13 nb^{-1}","p");
1008	<	leg2->AddEntry(hso["mc_qcd_clone"],"QCD","f");
1009	<	leg2->AddEntry(hso["mc_wj_clone"],"W+Jets","f");
1010	<	}
1011	<	else {
1012	<	leg2->AddEntry(hs["data_all_clone"],"Data   #int#font[12]{L}dt = 13 nb^{-1}","p");
1013	<	leg2->AddEntry(hs["mc_qcd_clone"],"QCD","f");
1014	<	leg2->AddEntry(hs["mc_wj_clone"],"W+Jets","f");
1007	>	//leg2->AddEntry(hso["data_all_clone"],"Data   #int#font[12]{L}dt = 35 nb^{-1}","p");
1008	>	leg2->AddEntry(hso["data_all_clone"],"Data","lp");
1009	>	//leg2->AddEntry(f2,TString("Extrapolation of "+mytxt),"l");
1010	>	leg2->AddEntry(f2,"Extrapolation","l");
1011	>	leg2->AddEntry(hso["mc_wj_clone"],style.WJetsText,"f");
1012	>	leg2->AddEntry(hso["mc_zj_clone"],style.DYZJetsText,"f");
1013	>	leg2->AddEntry(hso["mc_qcd_clone"],style.QCDText,"f");
1014	>	}
1015	>	else {
1016	>	//leg2->AddEntry(hs["data_all_clone"],"Data   #int#font[12]{L}dt = 35 nb^{-1}","p");
1017	>	leg2->AddEntry(hs["data_all_clone"],"Data","lp");
1018	>	//leg2->AddEntry(f2,TString("Extrapolation of "+mytxt),"l");
1019	>	leg2->AddEntry(f2,"Extrapolation","l");
1020	>	leg2->AddEntry(hs["mc_wj_clone"],style.WJetsText,"f");
1021	>	leg2->AddEntry(hs["mc_zj_clone"],style.DYZJetsText,"f");
1022	>	leg2->AddEntry(hs["mc_qcd_clone"],style.QCDText,"f");
1023		}
1024		leg2->Draw();
1025
#	Line 785 | Line 1027 | void fitRelIso_Data()
1027		//   fb2->SetParameters(par0[0],par0[1],par0[2]);
1028		//   fb2->SetLineColor(5);
1029		//   fb2->Draw("same");
1030	+	f2->Draw("same");
1031	+	fin->Draw("same");
1032	+	//hs["Data"]->GetXaxis()->SetRangeUser(0.,0.06);
1033	+	//hs["Data"]->GetYaxis()->SetRangeUser(0.3,300.);
1034	+	hs["Data"]->SetLineWidth(2);
1035	+	cvs[TString("cv1"+suffix2.str())]->SetLogy();
1036	+	hs["Data"]->DrawClone("esame"); // Draw on top of others
1037
1038		plotFile = plotFile(0,plotFile.Length()-4) + "_zoomed.pdf";
1039	<	//cvs[TString("cv1"+suffix2.str())]->Print(plotFile);
1039	>	//if (ibin == 2 || ibin == 3)
1040	>	cvs[TString("cv1"+suffix2.str())]->Print(plotFile);
1041
1042		hs.clear();
1043		hso.clear();
#	Line 800 | Line 1050 | void fitRelIso_Data()
1050		float average[2] = {0.,0.};
1051		float errors = 0.;
1052
1053	<	int nmax = sizeof(res_nbkg)/sizeof(int);
1054	<	for (int ii=0; ii<nmax; ii++) {
1055	<	sum[0] += res_nbkg[ii];
1056	<	sum[1] += res_nsig[ii];
1053	>	if ( nfiles > 1) {
1054	>	for (int ii=0; ii<nfiles; ii++) {
1055	>	sum[0] += res_nbkg[ii];
1056	>	sum[1] += res_nsig[ii];
1057	>	}
1058	>	average[0] = round(sum[0] / nfiles,0);
1059	>	average[1] = round(sum[1] / nfiles,0);
1060	>	errors = stdDev(res_nbkg,nfiles,average[0]);
1061	>
1062	>	cout << "Mean bkg estimated = " << average[0] << endl;
1063	>	cout << "Mean sig estimated = " << average[1] << endl;
1064	>	cout << "Standard deviation = " << errors << endl;
1065	>	cout << "Uncertainties      = " << errors/average[0]*100. << "% " << endl;
1066	>
1067	>	outTxtFile << "Estimated bkg      = ";
1068	>	for (int jj=0; jj<nfiles; ++jj)
1069	>	outTxtFile << res_nbkg[jj] << "  ";
1070	>	outTxtFile << endl;
1071	>	outTxtFile << "Mean bkg estimated = " << average[0] << endl;
1072	>	outTxtFile << "Mean sig estimated = " << average[1] << endl;
1073	>	outTxtFile << "Standard deviation = " << errors << endl;
1074	>	outTxtFile << "Uncertainties      = " << errors/average[0]*100. << "% " << endl;
1075		}
808	–	average[0] = round(sum[0] / nmax,0);
809	–	average[1] = round(sum[1] / nmax,0);
810	–	errors = stdDev(res_nbkg,nmax,average[0]);
811	–
812	–	cout << "Mean bkg estimated = " << average[0] << endl;
813	–	cout << "Mean sig estimated = " << average[1] << endl;
814	–	cout << "Standard deviation = " << errors << endl;
815	–	cout << "Uncertainties      = " << errors/average[0]*100. << "% " << endl;
816	–
817	–	outTxtFile << "Estimated bkg      = ";
818	–	for (int jj=0; jj<nmax;++jj)
819	–	outTxtFile << res_nbkg[jj] << "  ";
820	–	outTxtFile << endl;
821	–	outTxtFile << "Mean bkg estimated = " << average[0] << endl;
822	–	outTxtFile << "Mean sig estimated = " << average[1] << endl;
823	–	outTxtFile << "Standard deviation = " << errors << endl;
824	–	outTxtFile << "Uncertainties      = " << errors/average[0]*100. << "% " << endl;
825	–
1076		outTxtFile.close();
1077		}

Diff Legend

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