Annotation of ttbar/p20_taujets_note/tau_jettrig/trig_tau.C, revision 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 trig_tau.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