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();
}