ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing UserCode/Jeng/scripts/fitOldRelIso.C (file contents):
Revision 1.2 by jengbou, Mon Apr 27 22:38:45 2009 UTC vs.
Revision 1.4 by jengbou, Fri Oct 30 17:55:47 2009 UTC

#	Line 17 | Line 17 | using namespace std;
17		#include "Math/GSLIntegrator.h"
18		#include "Math/WrappedTF1.h"
19
20	<	// TString fileName = "skimmed/QCDbin4_OldRelIso_rebin.root";
21	<	TString fileName = "skimmed/QCDbin4_OldRelIso_FastSim.root";
20	>	// FullSim
21	>	// TString fileName = "skimmed226/QCDbin1_OldRelIso_rebin_TandL.root";
22	>	TString fileName = "skimmed226/QCDbin4_OldRelIso_wErrors_TandL_1.root";
23	>	// FastSim
24	>	// TString fileName = "skimmed226/QCDbin2_OldRelIso_FastSim_095.root";
25	>	// TString fileName = "skimmed226/QCDbin1_OldRelIso_wErrors_TandL_FastSim_0.root";
26
27	+	// Handles:
28	+	bool dynamic = true;
29		// fit region weight
30	<	double wxmin = 0.07; // 0.13 for FullSim 0.07 for FsatSim
31	<	double wxmax = 0.0;
30	>	double wxmax = 1.825; //2. FastSim ; 1.825 Full
31	>	double wxmin = wxmax/3.5; // 3.5
32		// Select fitting function:
33		TString s0 = "gaus";
34		// Old reliso:
35		TString isoname = "Old";
36	+	bool drawtail = false;
37	+	bool debug = false;
38
39		void fitOldRelIso()
40		{
41	<	gStyle->SetOptFit(01111);
42	<
41	>	gROOT->SetStyle("CMS");
42	>	//   gStyle->SetOptStat("n");
43	>	//   gStyle->SetOptTitle(1);
44	>	gStyle->SetOptFit(11111);
45	>
46		TFile *file = TFile::Open(fileName);
47		TTree *tree = (TTree*)file->Get("myTree");
48		TH1D *h0 = (TH1D*)file->Get("histos/h_RelIso_all");
49		TH1D *h1 = (TH1D*)file->Get("histos/h_RelIso_qcd");
50		Int_t jbin;
51	<	Double_t na;
51	>	Double_t na,na1;
52		tree->SetBranchAddress("jbin",&jbin);
53		tree->SetBranchAddress("na",&na);
54	+	if (fileName.Contains("TandL"))
55	+	tree->SetBranchAddress("na1",&na1);
56		tree->GetEntry(0);
57	<
57	>
58		// define fit region:
59		double ax0 = 0.0;
60		double axf = 0.0;
61		double bx0 = 0.0;
62		double bxf = 0.0;
63	+	double loosecut = 0.;
64
65		cout<<"Use "<<isoname<<" RelIso"<<endl;
66
67	<	// old reliso
68	<	ax0 = 0.95;
69	<	axf = 1.0;
70	<	bx0 = 0.0;
71	<	bxf = 0.8;
72	<
73	<	double normchi2 = 0.;
74	<	double ns, ns2;
75	<	ns = ns2 = 0.;
67	>	if ( isoname.Contains("Old") ) {
68	>	// old
69	>	loosecut = 0.9;
70	>	ax0 = 0.95;
71	>	axf = 1.0;
72	>	bx0 = 0.3; // Must be smaller than x at peak
73	>	bxf = 0.9; // Should not overlap signal region
74	>	}
75	>	else if ( isoname.Contains("New") ) {
76	>	// new
77	>	ax0 = 0.0;
78	>	axf = 1./19.; // ~0.053
79	>	loosecut = 1./9.; // ~0.11
80	>	bx0 = 0.2;
81	>	bxf = 2.;
82	>	}
83	>
84	>	int niter = 1;
85	>	double pass[3] = {0.,0.,0.};
86	>	double bound[2] = {0.,0.};
87	>	double ns[2] = {0.,0.};
88	>	double ns2 = 0.;
89
90	<	TCanvas *cv = new TCanvas("cv","cv",800,600);
90	>	TCanvas *cv = new TCanvas("cv","cv",800,800);
91
92		// Fit QCD in background region
93		cout<<"Fitting with "<<s0<<" function!"<<endl;
#	Line 68 | Line 95 | void fitOldRelIso()
95		TH1D *h_qcd = (TH1D*)h1->Clone();
96		h_all->SetLineWidth(2);
97		h_all->SetMinimum(0);
98	<	if( jbin == 4)
99	<	h_all->SetTitle(isoname+" RelIso distribution (#geq 4 jets)");
73	<	else if ( jbin < 4 ) {
74	<	TString title = isoname+" RelIso distribution (";
75	<	title += jbin;
76	<	title += " jets)";
77	<	h_all->SetTitle(title);
78	<	}
79	<	else{
80	<	cout << "Wrong bin number!" << endl;
81	<	}
82	<	h_all->Draw();
83	<
98	>
99	>	// Get num of bin at peak and bin width
100		Int_t nbxMax = h_all->GetXaxis()->FindBin(1.0)+1;
101		h_all->GetXaxis()->SetRange(1,h_all->GetXaxis()->FindBin(bxf));
102	<	Int_t nbMax=h_all->GetMaximumBin();
103	<	cout << "Bin number at peak= " << h_all->GetMaximumBin() << endl;
102	>	Int_t nbMax = h_all->GetMaximumBin();
103	>	cout << ">>>>> Bin number at peak= " << nbMax << endl;
104		h_all->GetXaxis()->SetRange(1,nbxMax);
105		TAxis *axis0 = h_all->GetXaxis();
106		double bwidth = axis0->GetBinWidth(nbMax);
107		double nbin = axis0->GetLast();
108		axf += bwidth;
109	<	cout << "bwidth = " << bwidth << endl;
94	<	bx0 = bwidth*nbMax - wxmin;
109	>	cout << ">>>>> bwidth = " << bwidth << endl;
110
111	<	cout << "\n>>>>> Start fitting control region: " << bx0 << " to " << bxf << endl;
111	>	TString title = "RelIso";
112	>	if( jbin == 4)
113	>	title += " (#geq 4 jets)";
114	>	else if ( jbin == 1 ) {
115	>	title += " (";
116	>	title += jbin;
117	>	title += " jet)";
118	>	}
119	>	else if ( jbin < 4 ) {
120	>	title += " (";
121	>	title += jbin;
122	>	title += " jets)";
123	>	h_RelIso_all->SetTitle(title);
124	>	}
125	>	else if ( jbin == 5 ) {
126	>	h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
127	>	}
128	>	else {
129	>	cout<<"Don't mess around!!! jbin should be less than 6!"<<endl;
130	>	return;
131	>	}
132	>	//   gStyle->SetErrorX(0);
133	>	h_all->GetXaxis()->SetTitle(title);
134	>	h_all->GetXaxis()->SetLabelFont(42);
135	>	h_all->GetXaxis()->SetLabelSize(0.03);
136	>	h_all->GetXaxis()->SetTitleFont(42);
137	>	h_all->GetXaxis()->SetTitleSize(0.035);
138	>	h_all->GetXaxis()->SetTitleOffset(1.2);
139
140	+	h_all->GetYaxis()->SetTitle("Events ( L_{Int}= 20 pb^{-1} )");
141	+	h_all->GetYaxis()->SetLabelFont(42);
142	+	h_all->GetYaxis()->SetLabelSize(0.03);
143	+	h_all->GetYaxis()->SetTitleFont(42);
144	+	h_all->GetYaxis()->SetTitleSize(0.035);
145	+	h_all->GetYaxis()->SetTitleOffset(1.6);
146	+	h_all->GetXaxis()->SetRangeUser(0.,1.2);
147	+	h_all->Draw("e");
148	+
149	+	cout << "\n>>>>> Iteration step: "<< niter << endl;
150	+	cout << "\n>>>>> Start fitting control region: " << bx0 << " to " << bxf << endl;
151	+
152		h_all->Fit(s0,"","",bx0,bxf);
153
154		// Very first fit function
#	Line 102 | Line 156 | void fitOldRelIso()
156		Int_t npar = f0->GetNpar();
157		Double_t chi2 = f0->GetChisquare();
158		Int_t ndof = f0->GetNDF();
159	<	normchi2 = chi2/ndof;
160	<	cout << "\n>>> Norm chi2 = " << normchi2 << endl;
159	>	pass[0] = chi2/ndof; // Very first fit norm chi2. Fixed range.
160	>	cout << ">>>> Norm chi2 = " << pass[0] << endl;
161	>	if (!dynamic)
162	>	pass[1] = pass[0];
163	>	else
164	>	pass[1] = 999.;
165
166		cout << " >>> Peak at " << f0->GetMaximumX(bx0,bxf) << endl;
167
168		double * par0 = new double [npar];
169		double * parErr0 = new double [npar];
170	+
171		for (Int_t i=0;i<npar;i++) {
172	+	printf("%s = %g +- %g\n",
173	+	f0->GetParName(i),
174	+	f0->GetParameter(i),
175	+	f0->GetParError(i)
176	+	);
177		par0[i] = f0->GetParameter(i);
178		parErr0[i] = f0->GetParError(i);
179	<	//       cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
179	>	}
180	>	bound[0] = bx0;
181	>	bound[1] = bxf;
182	>	if(dynamic) {
183	>	bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2);
184	>	bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2);
185		}
186
187	+	double delta = pass[1]-pass[2];
188	+
189	+	while (dynamic && niter <= 20 && abs(delta) > 0.00001) {
190	+	if (niter > 1)
191	+	pass[1] = pass[2];
192	+	if (delta != 0){
193	+	niter++;
194	+	cout << "\n>>>>> Iteration step: "<< niter << endl;
195	+	cout << ">> Previous step norm. chi2 = " << pass[1] << endl;
196	+	cout << ">>>> Starting fitting new range (bx0, bxf) = (" << bx0;
197	+	cout << ", " << bxf << ")" << endl;
198	+
199	+	h_all->Fit(s0,"0","",bx0,bxf);
200	+	TF1 *fi = (TF1*)h_all->GetFunction(s0);
201	+	pass[2] = fi->GetChisquare()/fi->GetNDF();
202	+	cout << ">> Current step norm. chi2  = " << pass[2] << endl;
203	+	delta = pass[1]-pass[2];
204	+	cout << ">> Delta = " << delta << endl;
205	+
206	+	if (delta > 0) {
207	+	bound[0] = bx0;
208	+	bound[1] = bxf;
209	+	printf("Function has %i parameters. Chisquare = %g\n",
210	+	npar,
211	+	fi->GetChisquare());
212	+	for (Int_t i=0;i<npar;i++) {
213	+	printf("%s = %g +- %g\n",
214	+	fi->GetParName(i),
215	+	fi->GetParameter(i),
216	+	fi->GetParError(i)
217	+	);
218	+	par0[i] = fi->GetParameter(i);
219	+	parErr0[i] = fi->GetParError(i);
220	+	//       cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
221	+	}
222	+	bx0 = fi->GetParameter(1) - wxmin*fi->GetParameter(2);
223	+	bxf = fi->GetParameter(1) + wxmax*fi->GetParameter(2);
224	+	cout << ">> Range for next iteration (bx0, bxf) = (" << bx0;
225	+	cout << ", " << bxf << ")" << endl;
226	+	}
227	+	else if (delta < 0) {
228	+	cout << " Use previous fit parameters for extrapolation fit" << endl;
229	+	delta = 0;
230	+	}
231	+	else {
232	+	//      bound[0] = (bound[0]<bx0)?bound[0]:bx0;
233	+	//      bound[1] = (bound[1]>bxf)?bound[1]:bxf;
234	+	cout << ">> Fit converges." << endl;
235	+	cout << ">> Best fitting region = (" << bound[0];
236	+	cout << ", " << bound[1] << ")" << endl;
237	+	}
238	+	if (niter == 2)
239	+	pass[0]=pass[2]; // First dynamic fit norm. chi2
240	+	}
241	+	}
242	+	if(debug){
243	+	cout << ">>> First dynamic fit norm. chi2   = " << pass[0] << endl;
244	+	cout << ">>> Min dynamic norm. chi2         = " << pass[1] << endl;
245	+	cout << ">>> Last fitting region (bx0, bxf) = (" << bx0 << ", " << bxf << ")"<< endl;
246	+	cout << ">>> Best fitting region (bx0, bxf) = (" << bound[0] << ", ";
247	+	cout << bound[1] << ")"<< endl;
248	+	}
249	+
250		// Get number of events within fitted region
251		TAxis *axis = h_all->GetXaxis();
252	<	int bmin = axis->FindBin(bx0);
253	<	int bmax = axis->FindBin(bxf);
252	>	int bmin = axis->FindBin(bound[0]);
253	>	int bmax = axis->FindBin(bound[1]);
254		double nfac1 = h_all->Integral(bmin,bmax);
255	<	nfac1 -= (h_all->GetBinContent(bmin))*(bx0-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
256	<	nfac1 -= (h_all->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bxf)/axis->GetBinWidth(bmax);
255	>	nfac1 -= (h_all->GetBinContent(bmin))*(bound[0]-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
256	>	nfac1 -= (h_all->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bound[1])/axis->GetBinWidth(bmax);
257		cout << ">>> Final Nfac = " << nfac1 << endl;
258
259		// ////////////////////////
260		// Histos and extrapolation
261	<
261	>	// ////////////////////////
262		TF1 *f1 = (TF1*)f0->Clone();
263	<	f1->SetRange(bx0,bxf);
263	>	f1->SetRange(bound[0],bound[1]);
264		for (int i=0;i<npar;i++) {
265		cout<<" >> Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
266		f1->SetParameter(i,par0[i]);
267		}
136	–
268		f1->SetLineColor(2);
269	+
270	+	h_qcd->SetLineWidth(2);
271		h_qcd->SetLineColor(40);
272		h_qcd->SetFillStyle(3018);
273		h_qcd->SetFillColor(38);
274	<	h_qcd->Draw("same");
274	>	h_qcd->Draw("same hist");
275		f1->Draw("same");
276
277		// Connect fitted and extrapolated region
145	–	cout << ">>> Extrapolation region: " << ax0 << " to " << axf << endl;
278		TF1 *fin = (TF1*)f1->Clone();
279	+	fin->SetRange(bound[1],ax0);
280		fin->SetLineStyle(2);
281		fin->SetLineColor(8);
149	–	fin->SetRange(bxf,ax0);
282		fin->Draw("same");
283
284	+	if (drawtail){
285	+	TF1 *fine = (TF1*)f1->Clone();
286	+	fine->SetRange(0.,bound[0]);
287	+	fine->SetLineStyle(2);
288	+	fine->SetLineColor(8);
289	+	fine->Draw("same");
290	+	}
291	+
292		TF1 *f2 = (TF1*)f1->Clone();
293		f2->SetRange(ax0,axf);
294		f2->SetLineColor(4);
#	Line 157 | Line 297 | void fitOldRelIso()
297		double *x=new double[np];
298		double *w=new double[np];
299		f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
300	<	ns = f2->IntegralFast(np,x,w,ax0,axf);
301	<	ns/=(f2->Integral(bx0,bxf));
302	<	ns*=nfac1;
303	<	cout<<">>> Ns by usual integral = "<<ns<<endl;
300	>
301	>	if (fileName.Contains("095") ||
302	>	fileName.Contains("TandL") ||
303	>	fileName.Contains("Tight")) {
304	>	ns[0] = f2->IntegralFast(np,x,w,ax0,axf);
305	>	ns[0]/=(f2->Integral(bound[0],bound[1]));
306	>	ns[0]*=nfac1;
307	>	}
308	>	else if (fileName.Contains("090") ||
309	>	fileName.Contains("Loose")) {
310	>	ns[0] = f2->IntegralFast(np,x,w,loosecut,axf);
311	>	ns[0]/=(f2->Integral(bound[0],bound[1]));
312	>	ns[0]*=nfac1;
313	>	}
314	>
315	>	ns[1] = f2->IntegralFast(np,x,w,loosecut,axf);
316	>	ns[1]/=(f2->Integral(bound[0],bound[1]));
317	>	ns[1]*=nfac1;
318	>
319	>	if (fileName.Contains("TandL") || fileName.Contains("090"))
320	>	cout << ">>> Loose Ns by usual integral    = " << ns[1] << endl;
321	>	else
322	>	cout << ">>> Tight Ns by usual integral    = " << ns[0] << endl;
323	>
324		delete [] x;
325		delete [] w;
326
#	Line 169 | Line 329 | void fitOldRelIso()
329		ROOT::Math::GSLIntegrator ig(1.E-8,1.E-8,1000);
330		ROOT::Math::WrappedTF1 wf(*g);
331		ig.SetFunction(wf);
332	<	ns2 = ig.Integral(ax0,axf);
333	<	ns2/=(g->Integral(bx0,bxf));
332	>	if (fileName.Contains("095") ||
333	>	fileName.Contains("TandL") ||
334	>	fileName.Contains("Tight"))
335	>	ns2 = ig.Integral(ax0,axf);
336	>	else if (fileName.Contains("090") ||
337	>	fileName.Contains("Loose"))
338	>	ns2 = ig.Integral(loosecut,axf);
339	>
340	>	ns2/=(g->Integral(bound[0],bound[1]));
341		ns2*=nfac1;
342	<	cout<<">>> Ns by MathMore integral = "<<ns2<<endl;
343	<
344	<	cv->Update();
345	<
346	<	TPaveStats *p1 = (TPaveStats*)h_all->GetListOfFunctions()->FindObject("stats");
347	<	p1->SetTextColor(kBlue);
348	<	p1->SetX1NDC(0.6);
342	>	cout<<">>> Ns by MathMore integral       = "<<ns2<<endl;
343	>
344	>	if(!dynamic){
345	>	cv->Update();
346	>	TPaveStats *p1 = (TPaveStats*)cv->GetPrimitive("stats");
347	>	}
348	>	else {
349	>	h_all->Fit(s0,"0","",bound[0],bound[1]);
350	>	cv->Update();
351	>	TPaveStats *p1 = (TPaveStats*)cv->GetPrimitive("stats");
352	>	}
353	>	p1->SetName("Gaussian fit");
354	>	TList *list = p1->GetListOfLines();
355	>	//   TText *tconst = p1->GetLineWith("Constant");
356	>	//   list->Remove(tconst);
357	>	TLatex *myt = new TLatex(0,0,"Gaussisn fit");
358	>	list->AddFirst(myt);
359	>	h_all->SetStats(0);
360	>	cv->Modified();
361	>
362	>	p1->SetTextFont(42);
363	>	p1->SetFillColor(0);
364	>	p1->SetFillStyle(0);
365	>	p1->SetBorderSize(0);
366	>	p1->SetX1NDC(0.65);
367		p1->SetX2NDC(0.88);
368	<	p1->SetY1NDC(0.62);
369	<	p1->SetY2NDC(0.88);
370	<	p1->Draw();
371	<
368	>	p1->SetY1NDC(0.65);
369	>	p1->SetY2NDC(0.85);
370	>	gPad->Update();
371	>
372		// Label histo
373	<	TLegend * leg1 = new TLegend(0.3,0.15,0.68,0.45);
373	>	TLegend * leg1 = new TLegend(0.22,0.65,0.47,0.85);
374	>	leg1->SetFillColor(0);
375	>	leg1->SetFillStyle(0);
376	>	leg1->AddEntry(f1,"Fit of all events in control region","l");
377	>	leg1->AddEntry(f2,"Extrapolation to signal region","l");
378		leg1->AddEntry(h_all,"All events (S+B)");
379		leg1->AddEntry(h_qcd,"QCD events");
191	–	leg1->AddEntry(f1,"Fit of all events in control region");
192	–	leg1->AddEntry(f2,"Extrapolation to signal region");
380		leg1->Draw();
381	+
382		// Print results
383		cout << "\n";
384		cout << ">>>>> Jet bin " << jbin << ":" << endl;
385	<	cout << ">>>>> Observed Ns = " << ns << endl;
385	>	cout << ">>>>> Bin number at peak = " << nbMax << endl;
386	>	cout << ">>>>> bwidth = " << bwidth << endl;
387	>	if (dynamic) {
388	>	cout << ">>>>> First dynamic fit norm. chi2   = " << pass[0] << endl;
389	>	cout << ">>>>> Min dynamic fit norm. chi2     = " << pass[1] << endl;
390	>	cout << ">>>>> Last fitting region (bx0, bxf) = (" << bx0 << ", " << bxf << ")"<< endl;
391	>	cout << ">>>>> Best fitting region (bx0, bxf) = (" << bound[0] << ", " << bound[1] << ")"<< endl;
392	>	}
393	>	cout << ">>>>> Observed Ns = " << ns[0] << endl;
394		cout << ">>>>> Expected Na = " << na << endl;
395	<	cout << "Deviation = " << (ns-na)/na*100.0 << " %" <<endl;
396	<
395	>	cout << "Deviation = " << (ns[0]-na)/na*100.0 << " %" <<endl;
396	>
397	>	if (fileName.Contains("TandL")) {
398	>	cout << "====================================" << endl;
399	>	cout << ">>>>> Observed Loose Ns = " << ns[1] << endl;
400	>	cout << ">>>>> Expected Loose Na = " << na1 << endl;
401	>	cout << "Deviation = " << (ns[1]-na1)/na1*100.0 << " %" <<endl;
402	>	}
403		}

Diff Legend

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