Annotation of ttbar/p20_taujets_note/tau_jettrig/tau_trf2D_fit.C, revision 1.1.1.1
1.1 uid12904 1: Double_t fit_ET(Double_t* xar, Double_t* par)
2: {
3: Double_t x = *xar;
4:
5: Double_t p2x = par[2]+x;
6: if (p2x==0) return 0;
7: //return par[0]*x + par[0]*exp(-par[1]*p2x);
8: //return par[0]*exp((x)/(p2x*p2x)) + par[1]*x;
9: return par[0]*exp((x)/(p2x*p2x)) + par[1]*x/(p2x*p2x);
10: //return par[0]*(x-par[1])/p2x ;
11: }
12:
13: Double_t fit_eta(Double_t* xar, Double_t* par)
14: {
15: Double_t x = *xar;
16:
17: Double_t dx = x - 1;
18: if (x<0) dx = x + 1;
19: //x = fabs(x);
20: dx = fabs(dx);
21: Double_t y = par[0] + par[1]*x + par[2]*x*x + par[3]*x*x*x + par[4]*x*x*x*x + par[5]*x*x*x*x*x + par[6]*x*x*x*x*x*x;
22: //if (dx<0.05) y = par[0]*exp(x) + par[1]*x*x + par[2]*x*x*x + par[3]*x*x*x*x;
23: return y;
24: }
25:
26: Double_t fit_eta_ET(Double_t* xar, Double_t* par)
27: {
28: Double_t ET = xar[0];
29: Double_t eta = xar[1];
30:
31: Double_t* par_ET = par+1;
32: Double_t* par_eta = par+4; // offset = 1 + npar_ET
33:
34: return par[0]*fit_ET(&ET,par_ET)*fit_eta(&eta,par_eta);
35: }
36:
37: Double_t sigma_ET(Double_t* xar, Double_t* par)
38: {
39: // given the convariance matrix of parameters, return error
40: Double_t x = *xar;
41:
42: Double_t p2x = par[2]+x;
43: if (p2x==0) return 1e10;
44:
45: Double_t dYdP0 = exp((x)/(p2x*p2x));
46: Double_t dYdP1 = x/(p2x*p2x);
47: Double_t dYdP2 = -2*par[0]*exp((x)/(p2x*p2x))*(x)/(p2x*p2x*p2x) - 2*par[1]*(x)/(p2x*p2x*p2x);
48:
49: // Double_t dYdP0 = exp((x)/p2x);
50: // Double_t dYdP1 = x;
51: // Double_t dYdP2 = -2*par[0]*exp((x)/p2x)*(x)/(p2x*p2x*p2x);
52:
53:
54: Double_t* p = par+3; // offset = npar_ET
55: Double_t s0 = *p++*dYdP0; s0 += *p++*dYdP1; s0 += *p++*dYdP2;
56: Double_t s1 = *p++*dYdP0; s1 += *p++*dYdP1; s1 += *p++*dYdP2;
57: Double_t s2 = *p++*dYdP0; s2 += *p++*dYdP1; s2 += *p++*dYdP2;
58:
59: Double_t sigma_y = dYdP0*s0 + dYdP1*s1 + dYdP2*s2;
60: //cout<<"dYdP0 :"<<dYdP0<<" dYdP1 :"<<dYdP1<<" dYdP2 :"<<dYdP2<<endl;
61: //cout<<"s0 :"<<s0<<" s1 :"<<s1<<" s2 :"<<s2<<endl;
62: //cout<<"sigma_y_ET :"<<sigma_y<<endl;
63: //sigma_y = 0.000005;
64: return sqrt(sigma_y);
65: }
66:
67: Double_t sigma_eta(Double_t* xar, Double_t* par)
68: {
69: // given the convariance matrix of parameters, return error
70: Double_t x = *xar;
71: //x = fabs(x);
72: //Double_t dYdP0 = 0;
73: //if (x!=0) dYdP0 = -par[0]/(x*x);
74: Double_t dYdP0 = 1;
75: Double_t dYdP1 = x*x;
76: Double_t dYdP2 = x*x*x;
77: Double_t dYdP3 = x*x*x*x;
78: Double_t dYdP4 = x*x*x*x*x;
79: Double_t dYdP5 = x*x*x*x*x*x;
80: Double_t dYdP6 = x*x*x*x*x*x*x;
81:
82:
83: Double_t* p = par+7; // offset = npar_eta
84: Double_t s0 = *p++*dYdP0; s0 += *p++*dYdP1; s0 += *p++*dYdP2; s0 += *p++*dYdP3; s0 += *p++*dYdP4; s0 += *p++*dYdP5; s0 += *p++*dYdP6;
85: Double_t s1 = *p++*dYdP0; s1 += *p++*dYdP1; s1 += *p++*dYdP2; s1 += *p++*dYdP3; s1 += *p++*dYdP4; s1 += *p++*dYdP5; s1 += *p++*dYdP6;
86: Double_t s2 = *p++*dYdP0; s2 += *p++*dYdP1; s2 += *p++*dYdP2; s2 += *p++*dYdP3; s2 += *p++*dYdP4; s2 += *p++*dYdP5; s2 += *p++*dYdP6;
87: Double_t s3 = *p++*dYdP0; s3 += *p++*dYdP1; s3 += *p++*dYdP2; s3 += *p++*dYdP3; s3 += *p++*dYdP4; s3 += *p++*dYdP5; s3 += *p++*dYdP6;
88: Double_t s4 = *p++*dYdP0; s4 += *p++*dYdP1; s4 += *p++*dYdP2; s4 += *p++*dYdP3; s4 += *p++*dYdP4; s4 += *p++*dYdP5; s4 += *p++*dYdP6;
89: Double_t s5 = *p++*dYdP0; s5 += *p++*dYdP1; s5 += *p++*dYdP2; s5 += *p++*dYdP3; s5 += *p++*dYdP4; s5 += *p++*dYdP5; s5 += *p++*dYdP6;
90: Double_t s6 = *p++*dYdP0; s6 += *p++*dYdP1; s6 += *p++*dYdP2; s6 += *p++*dYdP3; s6 += *p++*dYdP4; s6 += *p++*dYdP5; s6 += *p++*dYdP6;
91:
92: Double_t sigma_y = dYdP0*s0 + dYdP1*s1 + dYdP2*s2 + dYdP3*s3 + dYdP4*s4 + dYdP5*s5 + dYdP6*s6; //cout<<"sigma_y_eta :"<<sigma_y<<endl;
93: return sqrt(sigma_y);
94: }
95:
96: Double_t sigma_eta_ET(Double_t* xar, Double_t* par)
97: {
98: Double_t ET = xar[0];
99: Double_t eta = xar[1];
100:
101: Double_t* par_ET = par+1;
102: Double_t* par_eta = par+13; // offset = 1 + npar_ET + npar_ET*npar_ET
103:
104: Double_t val_ET = fit_ET(&ET,par_ET);
105: Double_t err_ET = sigma_ET(&ET,par_ET);
106: Double_t val_eta = fit_eta(&eta,par_eta);
107: Double_t err_eta = sigma_eta(&eta,par_eta);
108:
109: return par[0]*sqrt(val_ET*val_ET*err_eta*err_eta +
110: err_ET*err_ET*val_eta*val_eta);
111: }
112:
113: Double_t fit_ET_p(Double_t* xar, Double_t* par)
114: {
115: return fit_ET(xar,par) + sigma_ET(xar,par);
116: }
117:
118: Double_t fit_ET_m(Double_t* xar, Double_t* par)
119: {
120: return fit_ET(xar,par) - sigma_ET(xar,par);
121: }
122:
123: Double_t fit_eta_p(Double_t* xar, Double_t* par)
124: {
125: return fit_eta(xar,par) + sigma_eta(xar,par);
126: }
127:
128: Double_t fit_eta_m(Double_t* xar, Double_t* par)
129: {
130: return fit_eta(xar,par) - sigma_eta(xar,par);
131: }
132:
133: Char_t* next_name(Char_t* ch)
134: {
135: if (ch[3]!='9') { ch[3]+=1; return ch; }
136: ch[3] = '0';
137: if (ch[2]!='9') { ch[2]+=1; return ch; }
138: ch[2] = '0';
139: if (ch[1]!='9') { ch[1]+=1; return ch; }
140: ch[1] = '0';
141: return ch;
142: }
143:
144: void tau_trf2D_fit()
145: {
146: gROOT->SetStyle("Plain");
147: gStyle->SetOptStat(0);
148: gStyle->SetOptTitle(0);
149:
150: Int_t icol[4] = {1,2,4,3};
151:
152: Char_t ch_name[] = "t000";
153:
154: const Int_t nf = 1;
155:
156: TFile* parametr_file = new TFile("tautrf2D_highmetl_nobtag_2.root");
157:
158: Double_t hh_min[2] = {0.0, 0.0};
159: Double_t hh_max[2] = {1.05, 1.05};
160:
161: TH1F* hh_allj[2];
162: TH1F* hh_btag[2];
163: TH1F* hh_eff[2];
164: TF1* fitfun[2];
165: TF1* fitfun_p[2];
166: TF1* fitfun_m[2];
167:
168: TH2F* hh_allj_2d[nf];
169: TH2F* hh_btag_2d[nf];
170: TH2F* hh_allj_2d_tag[nf];
171: TH2F* hh_btag_2d_tag[nf];
172: TF2* fitfun_2d[nf];
173: TF2* fitfun_2d_err[nf];
174:
175:
176:
177: TH2F* TauTRFtau = new TH2F("TauJET_tau","Taus #eta vs pT", 18,20, 200, 40,-2, 2);
178: TH2F* TauTRFjet = new TH2F("TauJET_jet","Taus #eta vs pT", 18,20, 200, 40,-2, 2);
179:
180: // parameterize taggability on EMqcd, the nj=1 bin
181:
182: parametr_file->cd();
183: hh_allj[0] = (TH1F*)gROOT->FindObject("TauTRFjetPt")->Clone();
184: hh_allj[1] = (TH1F*)gROOT->FindObject("TauTRFjetEta")->Clone();
185: hh_btag[0] = (TH1F*)gROOT->FindObject("TauTRFtauPt")->Clone();
186: hh_btag[1] = (TH1F*)gROOT->FindObject("TauTRFtauEta")->Clone();
187:
188: hh_allj_2d[0] = (TH2F*)gROOT->FindObject("TauTRFjetEtaPt")->Clone();
189: hh_btag_2d[0] = (TH2F*)gROOT->FindObject("TauTRFtauEtaPt")->Clone();
190:
191: hh_allj_2d_tag[0] = (TH2F*)gROOT->FindObject("TauTRFjetEtaPtTAG")->Clone();
192: hh_btag_2d_tag[0] = (TH2F*)gROOT->FindObject("TauTRFtauEtaPtTAG")->Clone();
193:
194:
195: Char_t* xtit[2] = {"ET, GeV", "#eta"};
196:
197: const Int_t npar_ET = 3, npar_eta = 7;
198:
199: const Int_t npar_ET_pm = npar_ET + npar_ET*npar_ET;
200: const Int_t npar_eta_pm = npar_eta + npar_eta*npar_eta;
201:
202: TCanvas* c1 = new TCanvas("c1","c1",700,760);
203: c1->Divide(1,2);
204: for (Int_t k=0; k<2; ++k) {
205: c1->cd(k+1);
206: gPad->SetBottomMargin(0.15);
207: gPad->SetGridx();
208: gPad->SetGridy();
209:
210: Int_t n = hh_allj[k]->GetNbinsX();
211: Double_t xmin = hh_allj[k]->GetXaxis()->GetXmin();
212: Double_t xmax = hh_allj[k]->GetXaxis()->GetXmax();
213: hh_eff[k] = new TH1F(next_name(ch_name),"eff",n,xmin,xmax);
214: for (Int_t i=0; i<n; ++i) {
215: Double_t a = hh_allj[k]->GetBinContent(i+1);
216: Double_t b = hh_btag[k]->GetBinContent(i+1);
217: Double_t y = 0, e = 0;
218: if (a>0) {
219: y = b/a;
220: e = sqrt(b*(a-b)/a)/a;
221: }
222: hh_eff[k]->SetBinContent(i+1,y);
223: hh_eff[k]->SetBinError(i+1,e);
224: }
225: hh_eff[k]->SetMarkerStyle(20);
226: hh_eff[k]->SetMarkerSize(1.4);
227: hh_eff[k]->SetMarkerColor(2);
228: hh_eff[k]->SetMinimum(hh_min[k]);
229: hh_eff[k]->SetMaximum(hh_max[k]);
230: hh_eff[k]->GetXaxis()->SetTitle(xtit[k]);
231: hh_eff[k]->GetXaxis()->SetTitleSize(0.06);
232: hh_eff[k]->GetXaxis()->SetLabelSize(0.06);
233: hh_eff[k]->GetXaxis()->SetNdivisions(505);
234: hh_eff[k]->GetYaxis()->SetLabelSize(0.06);
235: hh_eff[k]->GetYaxis()->SetNdivisions(505);
236: hh_eff[k]->Draw("pe0");
237:
238: TMatrixD matrix_ET(npar_ET,npar_ET);
239: TMatrixD matrix_eta(npar_eta,npar_eta);
240: if (k==0) { // ET
241: fitfun[k] = new TF1(next_name(ch_name),fit_ET,xmin,xmax,npar_ET);
242: fitfun[k]->SetParameter(0,0.9);
243: fitfun[k]->SetParameter(1,0);
244: fitfun[k]->SetParameter(2,10);
245: hh_eff[k]->Fit(fitfun[k],"0");
246: cout<<"ET chi squared (per degree of freedom) :"<<fitfun[k]->GetChisquare()/24<<endl;
247: gMinuit->mnemat(matrix_ET.GetMatrixArray(),npar_ET);
248: } else { // eta
249: fitfun[k] = new TF1(next_name(ch_name),fit_eta,xmin,xmax,npar_eta);
250: fitfun[k]->SetParameter(0,0.9);
251: fitfun[k]->SetParameter(1,2e-03);
252: fitfun[k]->SetParameter(2,1.0e-03);
253: fitfun[k]->SetParameter(3,1.0e-03);
254: fitfun[k]->SetParameter(4,1.0e-03);
255: fitfun[k]->SetParameter(5,1.0e-03);
256: fitfun[k]->SetParameter(6,1.0e-03);
257: hh_eff[k]->Fit(fitfun[k],"0");
258: cout<<"Eta chi squared (per degree of freedom) :"<<fitfun[k]->GetChisquare()/46<<endl;
259: gMinuit->mnemat(matrix_eta.GetMatrixArray(),npar_eta);
260: }
261: fitfun[k]->SetLineColor(4);
262: fitfun[k]->SetLineWidth(1);
263: fitfun[k]->Draw("same");
264:
265: // draw error bands
266: Int_t np = 0;
267: TMatrixD* tm;
268: if (k==0) { // ET
269: np = npar_ET;
270: tm = &matrix_ET;
271: fitfun_p[k] = new TF1(next_name(ch_name),fit_ET_p,xmin,xmax,npar_ET_pm);
272: fitfun_m[k] = new TF1(next_name(ch_name),fit_ET_m,xmin,xmax,npar_ET_pm);
273: } else { // eta
274: np = npar_eta;
275: tm = &matrix_eta;
276: fitfun_p[k] = new TF1(next_name(ch_name),fit_eta_p,xmin,xmax,npar_eta_pm);
277: fitfun_m[k] = new TF1(next_name(ch_name),fit_eta_m,xmin,xmax,npar_eta_pm);
278: }
279: for (Int_t i = 0; i<np; ++i) {
280: fitfun_p[k]->SetParameter(i,fitfun[k]->GetParameter(i));
281: fitfun_m[k]->SetParameter(i,fitfun[k]->GetParameter(i));
282: }
283: Int_t ip = np;
284: for (Int_t i = 0; i<np; ++i) {
285: for (Int_t j = 0; j<np; ++j) {
286: Double_t covij = (*tm)(i,j);
287: fitfun_p[k]->SetParameter(ip,covij);
288: fitfun_m[k]->SetParameter(ip,covij);
289: ++ip;
290: }
291: }
292: fitfun_p[k]->SetLineColor(4);
293: fitfun_p[k]->SetLineWidth(1);
294: fitfun_p[k]->Draw("same");
295: fitfun_m[k]->SetLineColor(4);
296: fitfun_m[k]->SetLineWidth(1);
297: fitfun_m[k]->Draw("same");
298: }
299: c1->Modified();
300:
301: // 2d parameterization
302: TFile* ff_f2d = new TFile("f2d_tau.root","RECREATE");
303:
304: for (Int_t kf = 0; kf<nf; ++kf) {
305: Double_t sa = 0, sb = 0;
306: // fetch histogram parameters
307: Int_t nx = hh_allj_2d[kf]->GetNbinsX();
308: Double_t xmin = hh_allj_2d[kf]->GetXaxis()->GetXmin();
309: Double_t xmax = hh_allj_2d[kf]->GetXaxis()->GetXmax();
310: Int_t ny = hh_allj_2d[kf]->GetNbinsY();
311: Double_t ymin = hh_allj_2d[kf]->GetYaxis()->GetXmin();
312: Double_t ymax = hh_allj_2d[kf]->GetYaxis()->GetXmax();
313: // construct 2d function
314: Char_t f2d_name[256];
315: strcpy(f2d_name,"f2d_tau_");
316: //strcat(f2d_name,ss_name[kf]);
317: const Int_t npar = npar_ET+npar_eta+1;
318: Double_t pars[npar] = {1};
319: fitfun[0]->GetParameters(pars+1);
320: fitfun[1]->GetParameters(pars+npar_ET+1);
321: fitfun_2d[kf] = new TF2(f2d_name,fit_eta_ET,xmin,xmax,ymin,ymax,npar);
322: fitfun_2d[kf]->SetParameters(pars);
323: for (Int_t ix = 1; ix<=nx; ++ix) {
324: for (Int_t iy = 1; iy<=ny; ++iy) {
325: // Double_t a = hh_allj_2d[kf]->GetBinContent(ix,iy);
326: // Double_t b = hh_btag_2d[kf]->GetBinContent(ix,iy);
327: // Double_t xi = hh_allj_2d[kf]->GetXaxis()->GetBinCenter(ix);
328: // Double_t yi = hh_allj_2d[kf]->GetYaxis()->GetBinCenter(iy);
329: Double_t a = hh_allj_2d[kf]->GetBinContent(ix,iy);
330: Double_t b = hh_btag_2d[kf]->GetBinContent(ix,iy);
331: Double_t xi = hh_allj_2d[kf]->GetXaxis()->GetBinCenter(ix);
332: Double_t yi = hh_allj_2d[kf]->GetYaxis()->GetBinCenter(iy);
333: Double_t fi = fit_ET(&xi,pars+1)*fit_eta(&yi,pars+npar_ET+1);
334: TauTRFtau->Fill(xi,yi,fi);
335: if (a==0) TauTRFjet->Fill(xi,yi,0); else TauTRFjet->Fill(xi,yi,b/a);
336: //if (a!=0) cout<<"b/a :"<<b/a<<endl;
337: //cout<<"fi :"<<fi<<" a: "<<a<<" sa: "<<sa<<endl;
338: //cout<<"xi: "<<xi<<" yi: "<<yi<<endl;
339: sa += a*fi;
340: sb += b;
341: }
342: }
343: Double_t sf = 0;
344: cout << "sb[" << kf << "] = " << sb << endl;
345: cout << "sa[" << kf << "] = " << sa << endl;
346: if (sa>0) sf = (sb/sa);
347: //sf = 9.6;
348: cout << "sf[" << kf << "] = " << sf << endl;
349: fitfun_2d[kf]->SetParameter(0,sf);
350: fitfun_2d[kf]->Write();
351: // construct 2d function for error
352: strcat(f2d_name,"_err");
353: const Int_t npar1 = 1+npar_ET_pm+npar_eta_pm;
354: Double_t pars1[npar1];
355: pars1[0] = sf;
356: fitfun_p[0]->GetParameters(pars1+1);
357: fitfun_p[1]->GetParameters(pars1+npar_ET_pm+1);
358: fitfun_2d_err[kf] = new TF2(f2d_name,sigma_eta_ET,xmin,xmax,ymin,ymax,npar1);
359: fitfun_2d_err[kf]->SetParameters(pars1);
360: TauTRFjet->Write();TauTRFtau->Scale(sf); TauTRFtau->Write();
361: fitfun_2d_err[kf]->Write();
362: }
363: ff_f2d->Close();
364: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/small/back.gif){kind=icon}
Return to tau_trf2D_fit.C CVS log ![[TXT]](/icons/small/text.gif){kind=icon}
![[DIR]](/icons/small/dir.gif){kind=icon}
Up to [Nicadd] / ttbar / p20_taujets_note / tau_jettrig