#include "TROOT.h"
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TProfile.h"
#include "TF1.h"
#include "TGraph.h"

#include "interface/FuncDef.h"

void JetResol()
{

   setTDRStyle();
   gROOT->Reset();

   TFile *f1 = new TFile("JETResol.root","read");
   TTree *Analyze = (TTree*)f1->Get("Analysis");
  
   const Int_t numBin=7;
   double bins[numBin]={0,20,40,60,80,100,220};
   TProfile *ETSoverET = new TProfile("ETSoverET","Jet resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);

   double rms[numBin];
   double mean[numBin];

   TCanvas *c1 = new TCanvas("c1","JET resol",0,0,1000,650);
   c1->cd(); int frame=0;
   c1->Divide(6,2);

   float x[numBin-1];
   float y[numBin-1];
   float ex[numBin-1];
   float ey[numBin-1];

   float finval=0;//valeur a remplir puis a fitter

   for ( int i=0; i<(numBin-1); i++)  // première bin : i ==1 et pas i == 0
      {
        TAxis *xaxis = ETSoverET->GetXaxis();
        float binCenter = xaxis->GetBinCenter(i+1);
        int binMin = (int)xaxis->GetBinLowEdge(i+1);
        int binMax = (int)xaxis->GetBinUpEdge(i+1);
        cout<<"binMin "<<binMin<<" binMax "<<binMax<<" bin center "<<binCenter<<endl;
        char tempMin[500];
        if(i==0)binMin=5;
        sprintf(tempMin,"JetPTResol.NonSmearePT > %d",binMin);
        string mystringMin(tempMin);
        char tempMax[500];
        sprintf(tempMax,"JetPTResol.NonSmearePT < %d",binMax);
        string mystringMax(tempMax);
        char tempName[500];
        sprintf(tempName,"(JetPTResol.SmearePT)>>hETSoverET%d",i);
        string mystringName(tempName);

        c1->cd(++frame);
        GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
        x[i]=binCenter;
        finval=rms[i]/mean[i];
        y[i]=(finval*100);
        ex[i]=0;
        ey[i]=0;
        cout<<"finval "<<finval<<" mean val "<<mean[i]<<" rms value "<<rms[i]<<endl;
      }

   TCanvas *c2 = new TCanvas("c2","JET resol",100,100,600,450);
   c2->cd();

   TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,170);

   TGraphErrors *gr11 = new TGraphErrors((numBin-1),x,y,ex,ey);
   gr11->Draw("AP");
   gr11->SetTitle("");
   gr11->GetXaxis()->SetTitle("E_{T}^{MC} [GeV]");
   gr11->GetYaxis()->SetRangeUser(0,50);
   gr11->GetYaxis()->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{MC})_{fit}/<E_{T}^{rec}/E_{T}^{MC}>_{fit}");

   Double_t* params = fitfun->GetParameters();

   gr11->Fit("user","R");
   gr11->Fit("user","RI");
   gr11->Fit("user","RI");
   gr11->Fit("user","RI");
   char tempResol[500];
   sprintf(tempResol,"#frac{#sigma(E_{T}^{rec}/E_{T}^{MC})}{<E_{T}^{rec}/E_{T}^{MC}>} = #frac{%f}{#sqrt{E_{T}^{MC}}} #oplus #frac{%f}{E_{T}^{MC}} #oplus %f",params[1]/100,params[2]/100,params[0]/100);

   TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
   leg1->Draw();

   delete fitfun;
}


void ETmisResol()
{

   setTDRStyle();
   gROOT->Reset();

   TFile *f1 = new TFile("PTMIS.root","read");
   TTree *Analyze = (TTree*)f1->Get("Analysis");

   const Int_t numBin=6;
   double bins[numBin]={10,20,30,50,80,100};
   TProfile *ETSoverET = new TProfile("ETSoverET","ETmis resol",(numBin-1),bins,-1000,1000);

   double rms[numBin-1];
   double mean[numBin-1];

   TCanvas *c5 = new TCanvas("c5","gerrors2",0,0,1000,650);
   c5->cd(); int frame=0;
   c5->Divide(6,2);

   double x[numBin-1];
   double y[numBin-1];

   for ( int i=0; i<(numBin-1); i++)  // première bin : i ==1 et pas i == 0
      {
        TAxis *xaxis = ETSoverET->GetXaxis();
        float binCenter = xaxis->GetBinCenter(i+1);
        int binMin = (int)xaxis->GetBinLowEdge(i+1);
        int binMax = (int)xaxis->GetBinUpEdge(i+1);
        cout<<"binMin "<<binMin<<" binMax "<<binMax<<" bin center "<<binCenter<<endl;
        char tempMin[500];
        sprintf(tempMin,"ETmisResol.Et>%d",binMin);
        string mystringMin(tempMin);
        char tempMax[500];
        sprintf(tempMax,"ETmisResol.Et<%d",binMax);
        string mystringMax(tempMax);

        char tempName[500];
        sprintf(tempName,"ETmisResol.EtSmeare>>hETSoverET%d",i);
        string mystringName(tempName);

        c5->cd(++frame);
        GaussValuesETmis(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
        x[i]=binCenter;
        y[i]=rms[i];
        cout<<" mean val "<<mean[i]<<" rms value "<<rms[i]<<endl;

      }

   TCanvas *c6 = new TCanvas("c6","ETmis resolution",100,100,600,450);
   c6->cd();

   TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",0,200);

   TGraph *gr11 = new TGraph((numBin-1),x,rms);
   gr11->Draw("AP");
   gr11->GetXaxis()->SetTitle("MC^{MET} [GeV]");
   gr11->GetYaxis()->SetTitle("#sigma(REC^{MET}-MC^{MET})");
   gr11->Fit("user","RQ");
   gr11->Fit("user","RQ");
   gr11->Fit("user","RQ");
   gr11->Fit("user","R");

   TCanvas *c7 = new TCanvas("c7","ETmis resolution",100,100,600,450);
   c7->cd();

   TGraph *gr12 = new TGraph((numBin-1),x,mean);
   gr12->Draw("AP");
   gr12->GetXaxis()->SetTitle("MC^{MET} [GeV]");
   gr12->GetYaxis()->SetTitle("<(REC^{MET}-MC^{MET})>");


 delete fitfun;
}

void General()
{
  JetResol();
  ETmisResol();
}