Double_t fit_ET(Double_t* xar, Double_t* par)
{
Double_t x = *xar;
Double_t p2x = par[2]+x;
if (p2x==0) return 0;
//return par[0]*x + par[0]*exp(-par[1]*p2x);
//return par[0]*exp((x)/(p2x*p2x)) + par[1]*x;
return par[0]*exp((x)/(p2x*p2x)) + par[1]*x/(p2x*p2x);
//return par[0]*(x-par[1])/p2x ;
}
Double_t fit_eta(Double_t* xar, Double_t* par)
{
Double_t x = *xar;
Double_t dx = x - 1;
if (x<0) dx = x + 1;
//x = fabs(x);
dx = fabs(dx);
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;
//if (dx<0.05) y = par[0]*exp(x) + par[1]*x*x + par[2]*x*x*x + par[3]*x*x*x*x;
return y;
}
Double_t fit_eta_ET(Double_t* xar, Double_t* par)
{
Double_t ET = xar[0];
Double_t eta = xar[1];
Double_t* par_ET = par+1;
Double_t* par_eta = par+4; // offset = 1 + npar_ET
return par[0]*fit_ET(&ET,par_ET)*fit_eta(&eta,par_eta);
}
Double_t sigma_ET(Double_t* xar, Double_t* par)
{
// given the convariance matrix of parameters, return error
Double_t x = *xar;
Double_t p2x = par[2]+x;
if (p2x==0) return 1e10;
Double_t dYdP0 = exp((x)/(p2x*p2x));
Double_t dYdP1 = x/(p2x*p2x);
Double_t dYdP2 = -2*par[0]*exp((x)/(p2x*p2x))*(x)/(p2x*p2x*p2x) - 2*par[1]*(x)/(p2x*p2x*p2x);
// Double_t dYdP0 = exp((x)/p2x);
// Double_t dYdP1 = x;
// Double_t dYdP2 = -2*par[0]*exp((x)/p2x)*(x)/(p2x*p2x*p2x);
Double_t* p = par+3; // offset = npar_ET
Double_t s0 = *p++*dYdP0; s0 += *p++*dYdP1; s0 += *p++*dYdP2;
Double_t s1 = *p++*dYdP0; s1 += *p++*dYdP1; s1 += *p++*dYdP2;
Double_t s2 = *p++*dYdP0; s2 += *p++*dYdP1; s2 += *p++*dYdP2;
Double_t sigma_y = dYdP0*s0 + dYdP1*s1 + dYdP2*s2;
//cout<<"dYdP0 :"<<dYdP0<<" dYdP1 :"<<dYdP1<<" dYdP2 :"<<dYdP2<<endl;
//cout<<"s0 :"<<s0<<" s1 :"<<s1<<" s2 :"<<s2<<endl;
//cout<<"sigma_y_ET :"<<sigma_y<<endl;
//sigma_y = 0.000005;
return sqrt(sigma_y);
}
Double_t sigma_eta(Double_t* xar, Double_t* par)
{
// given the convariance matrix of parameters, return error
Double_t x = *xar;
//x = fabs(x);
//Double_t dYdP0 = 0;
//if (x!=0) dYdP0 = -par[0]/(x*x);
Double_t dYdP0 = 1;
Double_t dYdP1 = x*x;
Double_t dYdP2 = x*x*x;
Double_t dYdP3 = x*x*x*x;
Double_t dYdP4 = x*x*x*x*x;
Double_t dYdP5 = x*x*x*x*x*x;
Double_t dYdP6 = x*x*x*x*x*x*x;
Double_t* p = par+7; // offset = npar_eta
Double_t s0 = *p++*dYdP0; s0 += *p++*dYdP1; s0 += *p++*dYdP2; s0 += *p++*dYdP3; s0 += *p++*dYdP4; s0 += *p++*dYdP5; s0 += *p++*dYdP6;
Double_t s1 = *p++*dYdP0; s1 += *p++*dYdP1; s1 += *p++*dYdP2; s1 += *p++*dYdP3; s1 += *p++*dYdP4; s1 += *p++*dYdP5; s1 += *p++*dYdP6;
Double_t s2 = *p++*dYdP0; s2 += *p++*dYdP1; s2 += *p++*dYdP2; s2 += *p++*dYdP3; s2 += *p++*dYdP4; s2 += *p++*dYdP5; s2 += *p++*dYdP6;
Double_t s3 = *p++*dYdP0; s3 += *p++*dYdP1; s3 += *p++*dYdP2; s3 += *p++*dYdP3; s3 += *p++*dYdP4; s3 += *p++*dYdP5; s3 += *p++*dYdP6;
Double_t s4 = *p++*dYdP0; s4 += *p++*dYdP1; s4 += *p++*dYdP2; s4 += *p++*dYdP3; s4 += *p++*dYdP4; s4 += *p++*dYdP5; s4 += *p++*dYdP6;
Double_t s5 = *p++*dYdP0; s5 += *p++*dYdP1; s5 += *p++*dYdP2; s5 += *p++*dYdP3; s5 += *p++*dYdP4; s5 += *p++*dYdP5; s5 += *p++*dYdP6;
Double_t s6 = *p++*dYdP0; s6 += *p++*dYdP1; s6 += *p++*dYdP2; s6 += *p++*dYdP3; s6 += *p++*dYdP4; s6 += *p++*dYdP5; s6 += *p++*dYdP6;
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;
return sqrt(sigma_y);
}
Double_t sigma_eta_ET(Double_t* xar, Double_t* par)
{
Double_t ET = xar[0];
Double_t eta = xar[1];
Double_t* par_ET = par+1;
Double_t* par_eta = par+13; // offset = 1 + npar_ET + npar_ET*npar_ET
Double_t val_ET = fit_ET(&ET,par_ET);
Double_t err_ET = sigma_ET(&ET,par_ET);
Double_t val_eta = fit_eta(&eta,par_eta);
Double_t err_eta = sigma_eta(&eta,par_eta);
return par[0]*sqrt(val_ET*val_ET*err_eta*err_eta +
err_ET*err_ET*val_eta*val_eta);
}
Double_t fit_ET_p(Double_t* xar, Double_t* par)
{
return fit_ET(xar,par) + sigma_ET(xar,par);
}
Double_t fit_ET_m(Double_t* xar, Double_t* par)
{
return fit_ET(xar,par) - sigma_ET(xar,par);
}
Double_t fit_eta_p(Double_t* xar, Double_t* par)
{
return fit_eta(xar,par) + sigma_eta(xar,par);
}
Double_t fit_eta_m(Double_t* xar, Double_t* par)
{
return fit_eta(xar,par) - sigma_eta(xar,par);
}
Char_t* next_name(Char_t* ch)
{
if (ch[3]!='9') { ch[3]+=1; return ch; }
ch[3] = '0';
if (ch[2]!='9') { ch[2]+=1; return ch; }
ch[2] = '0';
if (ch[1]!='9') { ch[1]+=1; return ch; }
ch[1] = '0';
return ch;
}
void tau_trf2D_fit()
{
gROOT->SetStyle("Plain");
gStyle->SetOptStat(0);
gStyle->SetOptTitle(0);
Int_t icol[4] = {1,2,4,3};
Char_t ch_name[] = "t000";
const Int_t nf = 1;
TFile* parametr_file = new TFile("tautrf2D_highmetl_nobtag_2.root");
Double_t hh_min[2] = {0.0, 0.0};
Double_t hh_max[2] = {1.05, 1.05};
TH1F* hh_allj[2];
TH1F* hh_btag[2];
TH1F* hh_eff[2];
TF1* fitfun[2];
TF1* fitfun_p[2];
TF1* fitfun_m[2];
TH2F* hh_allj_2d[nf];
TH2F* hh_btag_2d[nf];
TH2F* hh_allj_2d_tag[nf];
TH2F* hh_btag_2d_tag[nf];
TF2* fitfun_2d[nf];
TF2* fitfun_2d_err[nf];
TH2F* TauTRFtau = new TH2F("TauJET_tau","Taus #eta vs pT", 18,20, 200, 40,-2, 2);
TH2F* TauTRFjet = new TH2F("TauJET_jet","Taus #eta vs pT", 18,20, 200, 40,-2, 2);
// parameterize taggability on EMqcd, the nj=1 bin
parametr_file->cd();
hh_allj[0] = (TH1F*)gROOT->FindObject("TauTRFjetPt")->Clone();
hh_allj[1] = (TH1F*)gROOT->FindObject("TauTRFjetEta")->Clone();
hh_btag[0] = (TH1F*)gROOT->FindObject("TauTRFtauPt")->Clone();
hh_btag[1] = (TH1F*)gROOT->FindObject("TauTRFtauEta")->Clone();
hh_allj_2d[0] = (TH2F*)gROOT->FindObject("TauTRFjetEtaPt")->Clone();
hh_btag_2d[0] = (TH2F*)gROOT->FindObject("TauTRFtauEtaPt")->Clone();
hh_allj_2d_tag[0] = (TH2F*)gROOT->FindObject("TauTRFjetEtaPtTAG")->Clone();
hh_btag_2d_tag[0] = (TH2F*)gROOT->FindObject("TauTRFtauEtaPtTAG")->Clone();
Char_t* xtit[2] = {"ET, GeV", "#eta"};
const Int_t npar_ET = 3, npar_eta = 7;
const Int_t npar_ET_pm = npar_ET + npar_ET*npar_ET;
const Int_t npar_eta_pm = npar_eta + npar_eta*npar_eta;
TCanvas* c1 = new TCanvas("c1","c1",700,760);
c1->Divide(1,2);
for (Int_t k=0; k<2; ++k) {
c1->cd(k+1);
gPad->SetBottomMargin(0.15);
gPad->SetGridx();
gPad->SetGridy();
Int_t n = hh_allj[k]->GetNbinsX();
Double_t xmin = hh_allj[k]->GetXaxis()->GetXmin();
Double_t xmax = hh_allj[k]->GetXaxis()->GetXmax();
hh_eff[k] = new TH1F(next_name(ch_name),"eff",n,xmin,xmax);
for (Int_t i=0; i<n; ++i) {
Double_t a = hh_allj[k]->GetBinContent(i+1);
Double_t b = hh_btag[k]->GetBinContent(i+1);
Double_t y = 0, e = 0;
if (a>0) {
y = b/a;
e = sqrt(b*(a-b)/a)/a;
}
hh_eff[k]->SetBinContent(i+1,y);
hh_eff[k]->SetBinError(i+1,e);
}
hh_eff[k]->SetMarkerStyle(20);
hh_eff[k]->SetMarkerSize(1.4);
hh_eff[k]->SetMarkerColor(2);
hh_eff[k]->SetMinimum(hh_min[k]);
hh_eff[k]->SetMaximum(hh_max[k]);
hh_eff[k]->GetXaxis()->SetTitle(xtit[k]);
hh_eff[k]->GetXaxis()->SetTitleSize(0.06);
hh_eff[k]->GetXaxis()->SetLabelSize(0.06);
hh_eff[k]->GetXaxis()->SetNdivisions(505);
hh_eff[k]->GetYaxis()->SetLabelSize(0.06);
hh_eff[k]->GetYaxis()->SetNdivisions(505);
hh_eff[k]->Draw("pe0");
TMatrixD matrix_ET(npar_ET,npar_ET);
TMatrixD matrix_eta(npar_eta,npar_eta);
if (k==0) { // ET
fitfun[k] = new TF1(next_name(ch_name),fit_ET,xmin,xmax,npar_ET);
fitfun[k]->SetParameter(0,0.9);
fitfun[k]->SetParameter(1,0);
fitfun[k]->SetParameter(2,10);
hh_eff[k]->Fit(fitfun[k],"0");
cout<<"ET chi squared (per degree of freedom) :"<<fitfun[k]->GetChisquare()/24<<endl;
gMinuit->mnemat(matrix_ET.GetMatrixArray(),npar_ET);
} else { // eta
fitfun[k] = new TF1(next_name(ch_name),fit_eta,xmin,xmax,npar_eta);
fitfun[k]->SetParameter(0,0.9);
fitfun[k]->SetParameter(1,2e-03);
fitfun[k]->SetParameter(2,1.0e-03);
fitfun[k]->SetParameter(3,1.0e-03);
fitfun[k]->SetParameter(4,1.0e-03);
fitfun[k]->SetParameter(5,1.0e-03);
fitfun[k]->SetParameter(6,1.0e-03);
hh_eff[k]->Fit(fitfun[k],"0");
cout<<"Eta chi squared (per degree of freedom) :"<<fitfun[k]->GetChisquare()/46<<endl;
gMinuit->mnemat(matrix_eta.GetMatrixArray(),npar_eta);
}
fitfun[k]->SetLineColor(4);
fitfun[k]->SetLineWidth(1);
fitfun[k]->Draw("same");
// draw error bands
Int_t np = 0;
TMatrixD* tm;
if (k==0) { // ET
np = npar_ET;
tm = &matrix_ET;
fitfun_p[k] = new TF1(next_name(ch_name),fit_ET_p,xmin,xmax,npar_ET_pm);
fitfun_m[k] = new TF1(next_name(ch_name),fit_ET_m,xmin,xmax,npar_ET_pm);
} else { // eta
np = npar_eta;
tm = &matrix_eta;
fitfun_p[k] = new TF1(next_name(ch_name),fit_eta_p,xmin,xmax,npar_eta_pm);
fitfun_m[k] = new TF1(next_name(ch_name),fit_eta_m,xmin,xmax,npar_eta_pm);
}
for (Int_t i = 0; i<np; ++i) {
fitfun_p[k]->SetParameter(i,fitfun[k]->GetParameter(i));
fitfun_m[k]->SetParameter(i,fitfun[k]->GetParameter(i));
}
Int_t ip = np;
for (Int_t i = 0; i<np; ++i) {
for (Int_t j = 0; j<np; ++j) {
Double_t covij = (*tm)(i,j);
fitfun_p[k]->SetParameter(ip,covij);
fitfun_m[k]->SetParameter(ip,covij);
++ip;
}
}
fitfun_p[k]->SetLineColor(4);
fitfun_p[k]->SetLineWidth(1);
fitfun_p[k]->Draw("same");
fitfun_m[k]->SetLineColor(4);
fitfun_m[k]->SetLineWidth(1);
fitfun_m[k]->Draw("same");
}
c1->Modified();
// 2d parameterization
TFile* ff_f2d = new TFile("f2d_tau.root","RECREATE");
for (Int_t kf = 0; kf<nf; ++kf) {
Double_t sa = 0, sb = 0;
// fetch histogram parameters
Int_t nx = hh_allj_2d[kf]->GetNbinsX();
Double_t xmin = hh_allj_2d[kf]->GetXaxis()->GetXmin();
Double_t xmax = hh_allj_2d[kf]->GetXaxis()->GetXmax();
Int_t ny = hh_allj_2d[kf]->GetNbinsY();
Double_t ymin = hh_allj_2d[kf]->GetYaxis()->GetXmin();
Double_t ymax = hh_allj_2d[kf]->GetYaxis()->GetXmax();
// construct 2d function
Char_t f2d_name[256];
strcpy(f2d_name,"f2d_tau_");
//strcat(f2d_name,ss_name[kf]);
const Int_t npar = npar_ET+npar_eta+1;
Double_t pars[npar] = {1};
fitfun[0]->GetParameters(pars+1);
fitfun[1]->GetParameters(pars+npar_ET+1);
fitfun_2d[kf] = new TF2(f2d_name,fit_eta_ET,xmin,xmax,ymin,ymax,npar);
fitfun_2d[kf]->SetParameters(pars);
for (Int_t ix = 1; ix<=nx; ++ix) {
for (Int_t iy = 1; iy<=ny; ++iy) {
// Double_t a = hh_allj_2d[kf]->GetBinContent(ix,iy);
// Double_t b = hh_btag_2d[kf]->GetBinContent(ix,iy);
// Double_t xi = hh_allj_2d[kf]->GetXaxis()->GetBinCenter(ix);
// Double_t yi = hh_allj_2d[kf]->GetYaxis()->GetBinCenter(iy);
Double_t a = hh_allj_2d[kf]->GetBinContent(ix,iy);
Double_t b = hh_btag_2d[kf]->GetBinContent(ix,iy);
Double_t xi = hh_allj_2d[kf]->GetXaxis()->GetBinCenter(ix);
Double_t yi = hh_allj_2d[kf]->GetYaxis()->GetBinCenter(iy);
Double_t fi = fit_ET(&xi,pars+1)*fit_eta(&yi,pars+npar_ET+1);
TauTRFtau->Fill(xi,yi,fi);
if (a==0) TauTRFjet->Fill(xi,yi,0); else TauTRFjet->Fill(xi,yi,b/a);
//if (a!=0) cout<<"b/a :"<<b/a<<endl;
//cout<<"fi :"<<fi<<" a: "<<a<<" sa: "<<sa<<endl;
//cout<<"xi: "<<xi<<" yi: "<<yi<<endl;
sa += a*fi;
sb += b;
}
}
Double_t sf = 0;
cout << "sb[" << kf << "] = " << sb << endl;
cout << "sa[" << kf << "] = " << sa << endl;
if (sa>0) sf = (sb/sa);
//sf = 9.6;
cout << "sf[" << kf << "] = " << sf << endl;
fitfun_2d[kf]->SetParameter(0,sf);
fitfun_2d[kf]->Write();
// construct 2d function for error
strcat(f2d_name,"_err");
const Int_t npar1 = 1+npar_ET_pm+npar_eta_pm;
Double_t pars1[npar1];
pars1[0] = sf;
fitfun_p[0]->GetParameters(pars1+1);
fitfun_p[1]->GetParameters(pars1+npar_ET_pm+1);
fitfun_2d_err[kf] = new TF2(f2d_name,sigma_eta_ET,xmin,xmax,ymin,ymax,npar1);
fitfun_2d_err[kf]->SetParameters(pars1);
TauTRFjet->Write();TauTRFtau->Scale(sf); TauTRFtau->Write();
fitfun_2d_err[kf]->Write();
}
ff_f2d->Close();
}
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>