ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing UserCode/Jeng/scripts/fitNewRelIso.C (file contents):
Revision 1.1 by jengbou, Mon Apr 27 22:03:02 2009 UTC vs.
Revision 1.3 by jengbou, Fri Oct 30 17:55:46 2009 UTC

#	Line 14 | Line 14
14		#include "TH1D.h"
15		#include "TTree.h"
16		using namespace RooFit;
17	<
18	<
17	>	// .493/.2114 j1
18	>	//.773/.3257  j2
19	>	// .863/.353 j3
20	>	//.808/.294 j4
21		void fitNewRelIso()
22		{
21	–	TString fileName = "skimmed/QCDbin2_005.root";
23		//   TString fileName = "skimmed/QCDbinall_NewRelIso.root";
24		//   TString fileName = "skimmed/QCDbin4_NewRelIso_CD30_rebin2.root";
25		//   TString fileName = "skimmed/QCDbin3_NewRelIso_CD.root";
26	<	//TString fileName = "skimmed/QCDbin4_rebin.root";
27	<	//TString fileName = "skimmed/QCDbin4_NewRelIso_FastSim.root";
26	>	//   TString fileName = "skimmed/QCDbin4_rebin.root";
27	>	//   TString fileName = "skimmed/QCDbin4_NewRelIso_FastSim.root";
28	>	//   TString fileName = "skimmed226/QCDbin4_005.root";
29	>	//   TString fileName = "skimmed226/QCDbin4_005.root";
30	>	TString fileName = "skimmed226/QCDbin123_NewRelIso_FastSim_005.root";
31
32		// Dynamic(true) or static(false) range
33		bool dynamic = true;
#	Line 31 | Line 35 | void fitNewRelIso()
35		double epslon = 1.e-5;
36		double epslon1 = 1.e-6;
37		// Dynamic fit region weight
38	<	double wxmin = 1.7; // 1.81 1.7
39	<	double wxmax = .94; // 0.9 .95
38	>	double itune1 = .25; // 0.15 for 1-3; 0.25 for 4
39	>	double itune2 = 1.65;  //1.65
40	>	double itune3 = 1.25; //1.1 full 1.25 fast to limit upper fit range
41	>	double wxmin = 0.; // 1.81 1.7  (1.85 for fastsim 1.825 for full
42	>	double wxmax = 0.; // 0.9 .95  (/1.65 fast /2. full)
43	>	// Select fitting function:
44	>	//   TString s0 = "landau";
45	>	// Select New or Old reliso:
46	>	//   TString isoname = "New";
47
48		Double_t ax0=0.;
49		Double_t axf=0.05;
50	<	Double_t axF=0.11;
51	<	Double_t bx0=0.2;
52	<	Double_t bxf=2.;
50	>	Double_t axF=0.1;
51	>	Double_t bx0=0.2; // 0.4 for 4  0.2 for 123
52	>	Double_t bxf=2.;  // 2.0 for 4  1.0 for 123
53		int niter = 1;
54		int nsiter = 0;
55		Double_t temp = 0.;
56		Double_t tempb = 0.;
57	<	Double_t pass[3] = {0.,0.,0.};
57	>	Double_t para[3] = {0.,0.,0.};
58		Double_t bound[3] = {10e6.,0.,0.};
59	+	int bmin;
60	+	int bmax;
61	+	int ndof;
62
63		gStyle->SetOptStat(1);
64		gStyle->SetPalette(1);
#	Line 59 | Line 73 | void fitNewRelIso()
73		tree->SetBranchAddress("jbin",&jbin);
74		tree->SetBranchAddress("na",&na);
75		tree->GetEntry(0);
76	<
76	>	//   if(jbin != 0)
77	>	//     itune1 *= jbin;
78	>
79		TH1D *h_all = (TH1D*)h0->Clone();
80		TH1D *h_qcd = (TH1D*)h1->Clone();
81		Int_t nbxMax = h_all->GetXaxis()->FindBin(2.0);
#	Line 72 | Line 88 | void fitNewRelIso()
88		cout << "=========> bwidth = " << bwidth << endl;
89
90		cout << "\n>>>>> Iteration step: "<< niter << endl;
91	<
91	>
92		RooRealVar x("x","x",0.,2.);
93		RooRealVar y("y","y",0.,2.);
94		RooDataHist *dh0 = new RooDataHist("dh0","dh0",x,Import(*h0));
#	Line 88 | Line 104 | void fitNewRelIso()
104		// ---------------------
105
106		// Declare variables x,mean,sigma with associated name, title, initial value and allowed range
107	<	// Gaussian
107	>	// Exponential parameters
108	>	RooRealVar lamda("lamda","decay constant of exponential",1,-100.,100.);
109	>	// Gaussian parameters
110		RooRealVar mean("mean","mean of gaussian",1.5,0.,2.);
111		RooRealVar sigma("sigma","width of gaussian",1,0.,100.);
112	<	// Landau
112	>	// Landau parameters
113		RooRealVar MPV("MPV","mean of landau",1,0,2);
114		RooRealVar Sig("Sig","sigma of landau",1,0,2);
115		// Build gaussian p.d.f in terms of x,mean and sigma
116	+	RooExponential expo("expo","exponential PDF",x,lamda);
117		RooGaussian gauss("gauss","gaussian PDF",x,mean,sigma);
118		RooLandau landau("landau","landau PDF",x,MPV,Sig);
119	+	// Build model
120	+	//   RooRealVar elfrac("elfrac","fraction of Exponential and Landau",0.05,0.,1.);
121	+	//   RooRealVar glfrac("glfrac","fraction of Gaussian and Landau",0.,0.,1.);
122	+	//   RooAddPdf model("model","g+l",RooArgList(gauss,landau),glfrac);
123
124		// Construct plot frame in 'x'
125		RooPlot* xframe = x.frame(Title("Landau p.d.f. with data"));
#	Line 105 | Line 128 | void fitNewRelIso()
128		// P l o t   m o d e l   a n d   c h a n g e   p a r a m e t e r   v a l u e s
129		// ---------------------------------------------------------------------------
130		// Change the value of Parameters
131	+	//   sigma.setVal(2.);
132	+	//   MPV.setVal(0.8412);
133	+	//   Sig.setVal(0.314);
134	+	// Plot gauss in frame (i.e. in x) and draw frame on canvas
135	+	//   gauss.plotOn(xframe);
136	+	//   landau.plotOn(xframe,LineColor(kRed));
137
109	–	// G e n e r a t e   e v e n t s
110	–	// -----------------------------
111	–
138		// Make a second plot frame in x and draw both the
139		// data and the p.d.f in the frame
140	<	dh0->plotOn(xframe);
141	<	dh1->plotOn(xframe,MarkerColor(9));
142	<	dhy0->plotOn(yframe);
143	<	dhy1->plotOn(yframe,MarkerColor(9));
140	>	dh0->plotOn(xframe,MarkerStyle(23),MarkerColor(2),LineColor(2),XErrorSize(0));
141	>	dh1->plotOn(xframe,MarkerStyle(22),MarkerColor(9),LineColor(9),XErrorSize(0));
142	>	dhy0->plotOn(yframe,MarkerStyle(23),MarkerColor(2),LineColor(2),XErrorSize(0));
143	>	dhy1->plotOn(yframe,MarkerStyle(22),MarkerColor(9),LineColor(9),XErrorSize(0));
144
145		// F i t   m o d e l   t o   d a t a
146		// -----------------------------
147
148		// Fit pdf to data
123	–	//   MPV.setRange(0.4,1.);
124	–	//   Sig.setRange(0.1,0.5);
125	–	//   Sig.setConstant(kFALSE);
149		RooFitResult*rest = landau.fitTo(*dh0,Range("ctrlregion"),PrintLevel(-1),Save());
150		rest->Print("v");
151		niter++;
152	<	// Print values of mean and sigma (that now reflect fitted values and errors)
153	<
154	<	Double_t xb0 = MPV.getVal()-wxmin*Sig.getVal();
155	<	Double_t xbf = MPV.getVal()+wxmax*Sig.getVal();
133	<
134	<	RooAbsReal* nQCDsigw = landau.createIntegral(x,Range("signalregion"));
135	<	RooAbsReal* nQCDbkgw = landau.createIntegral(x,Range("ctrlregion"));
152	>	RooAbsReal* nQCDsigw = landau.createIntegral(x,NormSet(x),Range("signalregion"));
153	>	RooAbsReal* nQCDbkgw = landau.createIntegral(x,NormSet(x),Range("ctrlregion"));
154	>	nQCDsigw->Print();
155	>	nQCDbkgw->Print();
156		// Draw all frames on a canvas
157	<	TCanvas* c = new TCanvas("fitNewRelIso","fitNewRelIso",1000,600);
158	<	c->Divide(2,1);
159	<	c->cd(1);
160	<
161	<	for (;niter <= np && nsiter < 2; niter++){
157	>	Double_t xb0 = bx0;
158	>	Double_t xbf = bxf;
159	>	wxmin = (MPV.getVal() - itune1)/Sig.getVal();
160	>	wxmax = wxmin/itune2;
161	>
162	>	TAxis *axis = h0->GetXaxis();
163	>	if (dynamic) {
164	>	TCanvas* c = new TCanvas("fitNewRelIso","fitNewRelIso",1000,600);
165	>	c->Divide(2,1);
166	>	c->cd(1);
167	>
168		xb0 = MPV.getVal()-wxmin*Sig.getVal();
169		xbf = MPV.getVal()+wxmax*Sig.getVal();
170	<	cout << "\n>>>>> Iteration step: "<< niter << endl;
171	<	cout << "=========> xb0, xbf = " << xb0 << ", " << xbf << endl;
172	<	x.setRange("ctrlregion",xb0,xbf);
173	<
174	<	// Construct a chi^2 of the data and the model,
175	<	// which is the input probability density scaled
176	<	// by the number of events in the dataset
177	<	RooChi2Var chi2("chi2","chi2",landau,*dh0,Range("ctrlregion"));
178	<	// Use RooMinuit interface to minimize chi^2
179	<	RooMinuit m(chi2);
180	<	m.migrad();
181	<	m.hesse();
182	<
183	<	if ((fabs(chi2.getVal() - bound[0]) <= epslon) ||
184	<	(fabs(chi2.getVal() - temp) <= epslon) ||
185	<	(fabs(xb0 - tempb) < = epslon1 )
186	<	){
187	<	nsiter++;
188	<	cout << "=========> Delta chi2 converges for " << nsiter << "time(s)" << endl;
189	<	}
190	<	if (chi2.getVal() < bound[0]){
191	<	bound[0] = chi2.getVal();
192	<	bound[1] = xb0;
193	<	bound[2] = xbf;
194	<	pass[1] = MPV.getVal();
195	<	pass[2] = Sig.getVal();
196	<	nQCDsigw = landau.createIntegral(x,NormSet(x),Range("signalregion"));
197	<	nQCDbkgw = landau.createIntegral(x,NormSet(x),Range("ctrlregion"));
198	<	cout << "=========> New Min chi2 = " << chi2.getVal() <<endl;
199	<	cout << "=========> Fitting range = " << bound[1] << " to " << bound[2] << endl;
200	<	landau.plotOn(xframe,NormRange("ctrlregion"),Range("ctrlregion"),LineColor(niter));
201	<	landau.paramOn(xframe,Layout(0.28,0.75,0.42));
202	<	landau.plotOn(xframe,NormRange("ctrlregion"),Range("signalregion"),LineColor(niter));
203	<	xframe->Draw();
204	<	nQCDsigw->Print();
205	<	nQCDbkgw->Print();
170	>	if (xb0 < 2.*axf )
171	>	xb0 = 2.*axf;
172	>	if (xbf > itune3)
173	>	xbf = itune3;
174	>
175	>	for (;niter <= np && nsiter < 2; niter++){
176	>	xb0 = MPV.getVal()-wxmin*Sig.getVal();
177	>	xbf = MPV.getVal()+wxmax*Sig.getVal();
178	>	if (xb0 < 2.*axf )
179	>	xb0 = 2.*axf;
180	>	if (xbf > itune3)
181	>	xbf = itune3;
182	>
183	>	bmin = axis->FindBin(xb0);
184	>	bmax = axis->FindBin(xbf);
185	>	// number of bins of non-empty entries (>5) - (num params. of Landau+1)
186	>	ndof = bmax-bmin+1-4;
187	>
188	>	cout << "\n>>>>> Iteration step: "<< niter << endl;
189	>	cout << "=========> xb0, xbf = " << xb0 << ", " << xbf << endl;
190	>	//     x.setRange("ctrlregion",xb0.getVal(),xbf.getVal());
191	>	x.setRange("ctrlregion",xb0,xbf);
192	>	x.setRange("intregion",0.,xbf);
193	>	x.setRange("centregion",axf,xb0);
194	>
195	>	// Construct a chi^2 of the data and the model,
196	>	// which is the input probability density scaled
197	>	// by the number of events in the dataset
198	>	//     RooChi2Var chi2("chi2","chi2",landau,*dh0,Range(xb0.getVal(),xbf.getVal()));
199	>	RooChi2Var chi2("chi2","chi2",landau,*dh0,Range("ctrlregion"));
200	>	// Use RooMinuit interface to minimize chi^2
201	>	RooMinuit m(chi2);
202	>	m.migrad();
203	>	m.hesse();
204	>
205	>	if ((fabs(chi2.getVal()/ndof - bound[0]) <= epslon) ||
206	>	(fabs(chi2.getVal()/ndof - temp) <= epslon) ||
207	>	(fabs(xb0 - tempb) < = epslon1 )
208	>	){
209	>	nsiter++;
210	>	cout << "=========> Conditions converge for " << nsiter << " time(s)" << endl;
211	>	}
212	>	if (chi2.getVal()/ndof < bound[0]){
213	>	bound[0] = chi2.getVal()/ndof;
214	>	bound[1] = xb0;
215	>	bound[2] = xbf;
216	>	para[1] = MPV.getVal();
217	>	para[2] = Sig.getVal();
218	>	nQCDsigw = landau.createIntegral(x,NormSet(x),Range("signalregion"));
219	>	nQCDbkgw = landau.createIntegral(x,NormSet(x),Range("ctrlregion"));
220	>	cout << "=========> New Min Norm. chi2 = " << chi2.getVal()/ndof <<endl;
221	>	cout << "=========> Fitting range = " << bound[1] << " to " << bound[2] << endl;
222	>	landau.plotOn(xframe,NormRange("intregion"),Range("ctrlregion"),LineColor(niter+1));
223	>	landau.paramOn(xframe,Layout(0.28,0.75,0.42));
224	>	landau.plotOn(xframe,NormRange("intregion"),Range("signalregion"),LineColor(niter+1));
225	>	landau.plotOn(xframe,NormRange("intregion"),Range("centregion"),
226	>	LineColor(niter+1),LineStyle(4));
227	>	xframe->Draw();
228	>	nQCDsigw->Print();
229	>	nQCDbkgw->Print();
230	>	}
231	>	temp = chi2.getVal()/ndof;
232	>	tempb = xb0;
233		}
234	<	temp = chi2.getVal();
235	<	tempb = xb0;
234	>	cout << "=========> Norm. chi2 of the last iteration = " << chi2.getVal()/ndof << endl;
235	>	cout << "=========> Fit range of the last iteration = " << xb0 << " to " << xbf << endl;
236	>	cout << "=========> MPV = " << MPV.getVal() << "; Sig = " << Sig.getVal() << endl;
237	>	cout << "================================================================" << endl;
238	>
239	>	// plot
240	>	Double_t mpv1 = para[1];
241	>	Double_t sig1 = para[2];
242	>	RooRealVar MPV1("MPV1","MPV1",mpv1);
243	>	RooRealVar Sig1("Sig1","Sig1",sig1);
244	>	RooLandau landau1("landau1","landau1",y,MPV1,Sig1);
245	>	cout << "=========> Minimum Norm. chi2 = " << bound[0] << endl;
246	>	cout << "=========> Final fit range = " << bound[1] << " to " << bound[2] << endl;
247	>	cout << "=========> Calculate nQCDsigw and nQCDbkgw: " << endl;
248	>	cout << "=========> Final MPV = para[1] = " << para[1] << endl;
249	>	cout << "=========> Final SIG = para[2] = " << para[2] << endl;
250	>
251	>	y.setRange("signalregion",ax0,axf);
252	>	y.setRange("centralregion",axf,bound[1]);
253	>	y.setRange("ctrlregion",bound[1],bound[2]);
254	>	nQCDsigw = landau1.createIntegral(y,NormSet(y),Range("signalregion"));
255	>	nQCDbkgw = landau1.createIntegral(y,NormSet(y),Range("ctrlregion"));
256	>	nQCDsigw->Print();
257	>	nQCDbkgw->Print();
258	>
259	>	landau1.plotOn(yframe,NormRange("ctrlregion"),Range("ctrlregion"));
260	>	landau1.plotOn(yframe,NormRange("ctrlregion"),Range("centralregion"),
261	>	LineColor(3),LineStyle(kDashed));
262	>	landau1.plotOn(yframe,NormRange("ctrlregion"),Range("signalregion"),LineColor(2));
263	>	c->cd(2);
264	>	yframe->Draw();
265	>
266		}
267	<	// plot
268	<	Double_t mpv1 = pass[1];
269	<	Double_t sig1 = pass[2];
270	<	RooRealVar MPV1("MPV1","MPV1",mpv1);
271	<	RooRealVar Sig1("Sig1","Sig1",sig1);
272	<	RooLandau landau1("landau1","landau1",y,MPV1,Sig1);
273	<
274	<	y.setRange("signalregion",ax0,axf);
275	<	y.setRange("centralregion",axf,bound[1]);
276	<	y.setRange("ctrlregion",bound[1],bound[2]);
277	<
278	<	landau1.plotOn(yframe,NormRange("ctrlregion"),Range("ctrlregion"));
279	<	landau1.plotOn(yframe,NormRange("ctrlregion"),Range("centralregion"),
280	<	LineColor(3),LineStyle(kDashed));
281	<	landau1.plotOn(yframe,NormRange("ctrlregion"),Range("signalregion"),LineColor(2));
199	<	c->cd(2);
200	<	yframe->Draw();
267	>	else { // Fixed range
268	>
269	>	TCanvas* c = new TCanvas("fitNewRelIso","fitNewRelIso",800,600);
270	>	para[1] = MPV.getVal();
271	>	para[2] = Sig.getVal();
272	>	bound[1] = bx0;
273	>	bound[2] = bxf;
274	>	bmin = axis->FindBin(bound[1]);
275	>	bmax = axis->FindBin(bound[2]);
276	>	// number of bins of non-empty entries (>5) - (num params. of Landau+1)
277	>	bound[0] = rest->minNll();
278	>	cout << "=========> Final fit range = " << bound[1] << " to " << bound[2] << endl;
279	>	// Print values of mean and sigma (that now reflect fitted values and errors)
280	>	cout << "=========> Final MPV = para[1] = " << para[1] << endl;
281	>	cout << "=========> Final SIG = para[2] = " << para[2] << endl;
282
283	+	x.setRange("signalregion",ax0,axf);
284	+	x.setRange("centralregion",axf,bound[1]);
285	+	x.setRange("ctrlregion",bound[1],bound[2]);
286	+
287	+	landau.plotOn(xframe,NormRange("ctrlregion"),Range("ctrlregion"));
288	+	landau.plotOn(xframe,NormRange("ctrlregion"),Range("centralregion"),
289	+	LineColor(3),LineStyle(kDashed));
290	+	landau.plotOn(xframe,NormRange("ctrlregion"),Range("signalregion"),LineColor(2));
291	+	xframe->Draw();
292	+
293	+	}
294	+	//   h0->SetLineColor(9);
295	+	//   h0->SetLineStyle(2);
296	+	//   h0->SetLineWidth(2);
297	+	h0->Draw("same");
298	+	//   h1->SetLineStyle(3);
299	+	//   h1->SetLineWidth(2);
300	+	h1->SetLineColor(9); //40
301	+	//   h1->SetFillStyle(3020); // 3018
302	+	//   h1->SetFillColor(41); //38
303	+	h1->Draw("same");
304
203	–	cout << "=========> chi2 of the last iteration = " << chi2.getVal() << endl;
204	–	cout << "=========> Fit range of the last iteration = " << xb0 << " to " << xbf << endl;
205	–	cout << "================================================================" << endl;
206	–	cout << "=========> Minimum chi2 = " << bound[0] << endl;
207	–	cout << "=========> Final fit range = " << bound[1] << " to " << bound[2] << endl;
208	–	cout << "=========> MPV = " << pass[1] << "; Sig = " << pass[2] << endl;
209	–
210	–	// Print values of mean and sigma (that now reflect fitted values and errors)
211	–	//   mean.Print();
212	–	//   sigma.Print();
213	–	//   MPV.Print();
214	–	//   Sig.Print();
215	–	cout << "=========> Final MPV = pass[1] = " << pass[1] << endl;
216	–	cout << "=========> Final MPV = pass[2] = " << pass[2] << endl;
217	–
305		// Get number of events within fitted region ==> cross check
306		TAxis *axis = h0->GetXaxis();
307	<	int bmin = axis->FindBin(xb0);
308	<	int bmax = axis->FindBin(xbf);
307	>	bmin = axis->FindBin(xb0);
308	>	bmax = axis->FindBin(xbf);
309	>	if ( !dynamic )
310	>	cout << "=========> Minimized FCN = " << bound[0] << endl;
311		double nfac = h0->Integral(bmin,bmax);
312		nfac -= (h0->GetBinContent(bmin))*(xb0-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
313		nfac -= (h0->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-xbf)/axis->GetBinWidth(bmax);
314		cout << "=========> Final Nfac = " << nfac << endl;
315	<
316	<	nQCDsigw->Print();
228	<	nQCDbkgw->Print();
315	>	//   nQCDsigw->Print();
316	>	//   nQCDbkgw->Print();
317		Double_t nQCDsig = nfac * nQCDsigw->getVal()/nQCDbkgw->getVal();
318		cout << "=========> Number of observed QCD evts = " << nQCDsig << endl;
319		// Print results
320		cout << "\n";
321		cout << ">>>>> Jet bin " << jbin << ":" << endl;
322	+	if ( dynamic ){
323	+	cout << ">>>>>  itune1 = " << itune1 << "; ";
324	+	cout << "itune2 = " << itune2 << endl;
325	+	cout << ">>>>>  wxmin = " << wxmin << "; ";
326	+	cout << "wxmax = " << wxmax << endl;
327	+	}
328		cout << ">>>>> Observed Ns = " << nQCDsig << endl;
329		cout << ">>>>> Expected Na = " << na << endl;
330		cout << "Deviation = " << (nQCDsig-na)/na*100.0 << " %" <<endl;
331	<
331	>
332		}

Diff Legend

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