#include "MyAnalysis.h"
#include "CutAndCount.h"
#include "METSmearing.h"
#include "MakeTree.h"
#include "Plotter.h"
#include <iostream>
#include <TChain.h>
#include <TLegend.h>
#include <TH1.h>
#include <TGraphAsymmErrors.h>
#include <string>




/////////////////////////////////////////
// Functions to get the CMS paper plots 
/////////////////////////////////////////
// getPaperPlot: returns TH1 with CMS experimental results
// getPaperPlotTheo: returns TGraphAsymErrors with theory expeactations
// DrawPaperPlot: draws both plots

// Function to get the CMS paper plots (experimental results)
//=============================================================

TH1F* getPaperPlot(float ptCut){

  float binsPlot[9];
  binsPlot[0]=-0.1;
  binsPlot[1]=0.0;
  binsPlot[2]=0.4;
  binsPlot[3]=0.8;
  binsPlot[4]=1.2;
  binsPlot[5]=1.5;
  binsPlot[6]=1.8;
  binsPlot[7]=2.1;
  binsPlot[8]=2.7;

  float valuesFromPaper[6];
  float errorsFromPaper[6];

  if(ptCut==25){
    valuesFromPaper[0]=0.147;
    valuesFromPaper[1]=0.159;
    valuesFromPaper[2]=0.184;
    valuesFromPaper[3]=0.207;
    valuesFromPaper[4]=0.231;
    valuesFromPaper[5]=0.253;
    
    errorsFromPaper[0]=0.010;
    errorsFromPaper[1]=0.0092;
    errorsFromPaper[2]=0.0125;
    errorsFromPaper[3]=0.0122;
    errorsFromPaper[4]=0.0136;
    errorsFromPaper[5]=0.0161;
    
  }else{
    
    valuesFromPaper[0]=0.131;
    valuesFromPaper[1]=0.139;
    valuesFromPaper[2]=0.158;
    valuesFromPaper[3]=0.185;
    valuesFromPaper[4]=0.202;
    valuesFromPaper[5]=0.231;
    
    errorsFromPaper[0]=0.0122;
    errorsFromPaper[1]=0.0114;
    errorsFromPaper[2]=0.0147;
    errorsFromPaper[3]=0.0136;
    errorsFromPaper[4]=0.0144;
    errorsFromPaper[5]=0.0174;
  }
  
  TString nameh="AsymPap_";
  nameh+=ptCut;


  TH1F *histAsym = new TH1F(nameh.Data(),nameh.Data(),8,binsPlot);
  
  for (int i=0;i<6;i++){
    histAsym->SetBinContent(i+2,valuesFromPaper[i]);
    histAsym->SetBinError(i+2,errorsFromPaper[i]);
  }
  
  histAsym->SetTitle("");
  histAsym->SetXTitle("Lepton #eta");
  histAsym->SetYTitle("Lepton charge asymetry");  
  histAsym->SetMarkerStyle(21);
  histAsym->GetYaxis()->SetRangeUser(0.02,0.34);
  histAsym->GetYaxis()->CenterTitle();
  histAsym->GetXaxis()->CenterTitle();
  return histAsym;
  
}

// Function to get the CMS paper plots (theory expectations)
//=============================================================

vector<TGraphAsymmErrors*> getPaperPlotTheo(float ptCut){

  float binsPlot[7];
  binsPlot[0]=0.0;
  binsPlot[1]=0.4;
  binsPlot[2]=0.8;
  binsPlot[3]=1.2;
  binsPlot[4]=1.5;
  binsPlot[5]=1.8;
  binsPlot[6]=2.1;

  float x[6];
  float ex[6];
  for (int i=0;i<6;i++){
    x[i]=binsPlot[i]+0.5*(binsPlot[i+1]-binsPlot[i]);
    ex[i]=0;
  }


  float valuesFromPaper[2][6]; // 0 = RESBOS, 1 = MCFM
  float errorsLowFromPaper[2][6]; // 0 = RESBOS, 1 = MCFM
  float errorsHighFromPaper[2][6]; // 0 = RESBOS, 1 = MCFM

  if(ptCut==25){

    valuesFromPaper[0][0]=0.157;
    valuesFromPaper[0][1]=0.169;
    valuesFromPaper[0][2]=0.193;
    valuesFromPaper[0][3]=0.22;
    valuesFromPaper[0][4]=0.246;
    valuesFromPaper[0][5]=0.265;
    
    errorsHighFromPaper[0][0]=0.8;
    errorsHighFromPaper[0][1]=0.9;
    errorsHighFromPaper[0][2]=0.8;
    errorsHighFromPaper[0][3]=0.8;
    errorsHighFromPaper[0][4]=0.8;
    errorsHighFromPaper[0][5]=0.8;

    errorsLowFromPaper[0][0]=1.0;
    errorsLowFromPaper[0][1]=1.0;
    errorsLowFromPaper[0][2]=1.1;
    errorsLowFromPaper[0][3]=1.1;
    errorsLowFromPaper[0][4]=1.1;
    errorsLowFromPaper[0][5]=1.0;

    valuesFromPaper[1][0]=0.153;
    valuesFromPaper[1][1]=0.167;
    valuesFromPaper[1][2]=0.192;
    valuesFromPaper[1][3]=0.220;
    valuesFromPaper[1][4]=0.245;
    valuesFromPaper[1][5]=0.263;
    
    errorsHighFromPaper[1][0]=0.8;
    errorsHighFromPaper[1][1]=0.9;
    errorsHighFromPaper[1][2]=0.8;
    errorsHighFromPaper[1][3]=0.8;
    errorsHighFromPaper[1][4]=0.8;
    errorsHighFromPaper[1][5]=0.8;

    errorsLowFromPaper[1][0]=1.0;
    errorsLowFromPaper[1][1]=1.0;
    errorsLowFromPaper[1][2]=1.1;
    errorsLowFromPaper[1][3]=1.1;
    errorsLowFromPaper[1][4]=1.1;
    errorsLowFromPaper[1][5]=1.0;
    
  }else{
    
    valuesFromPaper[0][0]=0.134;
    valuesFromPaper[0][1]=0.146;
    valuesFromPaper[0][2]=0.169;
    valuesFromPaper[0][3]=0.196;
    valuesFromPaper[0][4]=0.222;
    valuesFromPaper[0][5]=0.245;
    
    errorsHighFromPaper[0][0]=0.7;
    errorsHighFromPaper[0][1]=0.8;
    errorsHighFromPaper[0][2]=0.8;
    errorsHighFromPaper[0][3]=0.8;
    errorsHighFromPaper[0][4]=0.8;
    errorsHighFromPaper[0][5]=0.8;

    errorsLowFromPaper[0][0]=0.9;
    errorsLowFromPaper[0][1]=0.8;
    errorsLowFromPaper[0][2]=1.0;
    errorsLowFromPaper[0][3]=1.0;
    errorsLowFromPaper[0][4]=1.1;
    errorsLowFromPaper[0][5]=1.1;

    valuesFromPaper[1][0]=0.131;
    valuesFromPaper[1][1]=0.145;
    valuesFromPaper[1][2]=0.168;
    valuesFromPaper[1][3]=0.194;
    valuesFromPaper[1][4]=0.219;
    valuesFromPaper[1][5]=0.241;
    
    errorsHighFromPaper[1][0]=0.7;
    errorsHighFromPaper[1][1]=0.8;
    errorsHighFromPaper[1][2]=0.8;
    errorsHighFromPaper[1][3]=0.8;
    errorsHighFromPaper[1][4]=0.8;
    errorsHighFromPaper[1][5]=0.8;

    errorsLowFromPaper[1][0]=0.9;
    errorsLowFromPaper[1][1]=0.8;
    errorsLowFromPaper[1][2]=1.0;
    errorsLowFromPaper[1][3]=1.0;
    errorsLowFromPaper[1][4]=1.1;
    errorsLowFromPaper[1][5]=1.1;

  }
  float yR[6];
  float yM[6];
  float eyRL[6];
  float eyML[6];
  float eyRH[6];
  float eyMH[6];
  
  for (int i=0;i<6;i++){
    yR[i]= valuesFromPaper[0][i];
    yM[i]= valuesFromPaper[1][i];
    eyRL[i]=errorsLowFromPaper[0][i]/100.;
    eyRH[i]=errorsHighFromPaper[0][i]/100.;
    eyML[i]=errorsLowFromPaper[1][i]/100.;
    eyMH[i]=errorsHighFromPaper[1][i]/100.;
  }

  TGraphAsymmErrors *histR = new TGraphAsymmErrors(6,x,yR,ex,ex,eyRL,eyRH);
  TGraphAsymmErrors *histM = new TGraphAsymmErrors(6,x,yM,ex,ex,eyML,eyMH);
  histR->SetFillColor(5);
  histM->SetFillColor(8);
  histR->SetLineColor(kBlue);
  histM->SetLineColor(kBlue);
  histR->SetLineStyle(kDashed);
  histM->SetLineStyle(kDashed);
  histR->SetLineWidth(2);
  histM->SetLineWidth(2);


  vector<TGraphAsymmErrors*> res;
  res.push_back(histR);
  res.push_back(histM);
  return res;

}

// Function to draw the CMS paper plots (exp+theory)
//==================================================

void SetStyle(){

  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);

  gStyle->SetTitleOffset(1.0,"x");
  gStyle->SetTitleOffset(1.2,"y");

  gStyle->SetPadBottomMargin(0.1);
  gStyle->SetPadTopMargin(0.05);
  gStyle->SetPadLeftMargin(0.10);
  gStyle->SetPadRightMargin(0.05);
  
}

void DrawPaperPlot(){

  SetStyle();
  
  TH1F *Paper_30  =  getPaperPlot(30);
  TH1F *Paper_25  =  getPaperPlot(25);

  vector<TGraphAsymmErrors*> theo25=getPaperPlotTheo(25);
  vector<TGraphAsymmErrors*> theo30=getPaperPlotTheo(30);

  Paper_25->SetMarkerColor(kBlack);
  Paper_25->SetMarkerColor(kBlack);
  Paper_30->SetMarkerColor(kBlack);
  Paper_30->SetMarkerColor(kBlack);
  Paper_25->SetLineColor(kBlack);
  Paper_30->SetLineColor(kBlack);

  TLegend *leg25= new TLegend(0.25,0.7,0.5,0.85); 
  leg25->SetBorderSize(0);
  leg25->SetFillStyle(0);
  leg25->SetTextFont(42);
  leg25->SetHeader("pT > 25 GeV");
  leg25->AddEntry(Paper_25,"CMS paper","lm"); 
  leg25->AddEntry(theo25[0],"Resbos","fl"); 
  leg25->AddEntry(theo25[1],"MCFM","fl"); 
  
  TLegend *leg30= new TLegend(0.25,0.7,0.5,0.85); 
  leg30->SetBorderSize(0);
  leg30->SetFillStyle(0);
  leg30->SetTextFont(42);
  leg30->SetHeader("pT > 30 GeV");
  leg30->AddEntry(Paper_30,"CMS paper","lm"); 
  leg30->AddEntry(theo30[0],"Resbos","fl"); 
  leg30->AddEntry(theo30[1],"MCFM","fl"); 
  
  TCanvas *can25 = new TCanvas("can25","can25");
  Paper_25->Draw("E1");
  theo25[0]->Draw("3same");
  theo25[0]->Draw("Lsame");
  theo25[1]->Draw("3same");
  theo25[1]->Draw("Lsame");
  Paper_25->Draw("E1Same");
  leg25->Draw("same");
  can25->Print("plots/Paper25.pdf","pdf");
  

  TCanvas *can30 = new TCanvas("can30","can30");
  Paper_30->Draw("E1");
  theo30[0]->Draw("3same");
  theo30[0]->Draw("Lsame");
  theo30[1]->Draw("3same");
  theo30[1]->Draw("Lsame");
  Paper_30->Draw("E1same");
  leg30->Draw("same");
  can30->Print("plots/Paper30.pdf","pdf");

}

////////////////////////////////
// Exercice 1.1: Do basic plots
////////////////////////////////

void BasicStackPlots(){


  float lumi = 50.;
   
   MyAnalysis *A = new MyAnalysis();
   TChain* ch = new TChain("events");
   ch->Add("files/data.root"); 
   ch->Process(A);
   
   MyAnalysis *B = new MyAnalysis();
   TChain* ch2 = new TChain("events");
   ch2->Add("files/ttbar.root");
   ch2->Process(B);
   
   MyAnalysis *C = new MyAnalysis();
   TChain* ch3 = new TChain("events");
   ch3->Add("files/wjets.root");
   ch3->Process(C);
   
   MyAnalysis *D = new MyAnalysis();
   TChain* ch4 = new TChain("events");
   ch4->Add("files/dy.root");
   ch4->Process(D);
   
   MyAnalysis *E = new MyAnalysis();
   TChain* ch5 = new TChain("events");
   ch5->Add("files/ww.root");
   ch5->Process(E);

   MyAnalysis *F = new MyAnalysis();
   TChain* ch6 = new TChain("events");
   ch6->Add("files/wz.root");
   ch6->Process(F);

   MyAnalysis *G = new MyAnalysis();
   TChain* ch7 = new TChain("events");
   ch7->Add("files/zz.root");
   ch7->Process(G);

   MyAnalysis *H = new MyAnalysis();
   TChain* ch8 = new TChain("events");
   ch8->Add("files/qcd.root");
   ch8->Process(H);
   
   MyAnalysis *I = new MyAnalysis();
   TChain* ch9 = new TChain("events");
   ch9->Add("files/single_top.root");
   ch9->Process(I);
   
   Plotter P;
   P.SetData(A->histograms, std::string("Data"));
   P.AddBg(B->histograms, std::string("TTbar"));
   P.AddBg(C->histograms, std::string("Wjets"));
   P.AddBg(D->histograms, std::string("DY"));
   P.AddBg(H->histograms, std::string("QCD"));
   P.AddBg(E->histograms, std::string("WW"));
   P.AddBg(F->histograms, std::string("WZ"));
   P.AddBg(G->histograms, std::string("ZZ"));
   P.AddBg(I->histograms, std::string("Single Top"));
  
   P.Plot(string("plots/BasicStackPlots.pdf"),false);
   P.Plot(string("plots/BasicStackPlotsLog.pdf"),true);

 
}

///////////////////////////////////////////////////
// Exercice 1.2: Compute raw cut and count asymetry
///////////////////////////////////////////////////

// Function getAsym: computes the raw cut and count asymetry for a given MetCut and PtCut using class CutAndCount
// ================================================================================================================

TH1F* getAsym(float MetCut, float PtCut, bool doW){

  // Compute asymetry for pt cut 25
  CutAndCount *A = new CutAndCount();
  TChain* ch = new TChain("events");
  ch->Add("files/data.root"); 
  A->SetPtCut(PtCut);
  A->SetMetCut(MetCut);
  if(!doW)  A->SetMetCut(0);
  A->SetDoW(doW);
  ch->Process(A);

  vector<float> etaLow=A->EtaLow;
  vector<float> etaHigh=A->EtaHigh;
  
  vector<float> asym=A->Asymetry;
  vector<float> error=A->AsymetryError;
  
 
  int nbins=6;
  float bins[7];
  float binsPlot[9];
  binsPlot[0]=-0.1;
  binsPlot[8]=2.7;

  for (int i=0;i<etaLow.size();i++){
    cout<<" Asym["<<etaLow[i]<<"<eta<"<<etaHigh[i]<<"] = "<<asym[i]<< " +/- "<< error[i]<< endl;
    bins[i]=etaLow[i];
    binsPlot[i+1]=bins[i];
  }
  bins[nbins]=etaHigh[nbins-1];
  binsPlot[nbins+1]=bins[nbins];
 
  TString nameh="Asym_";
  nameh+=PtCut;
  nameh+="_";
  nameh+=MetCut;

  TH1F *histAsym = new TH1F(nameh.Data(),nameh.Data(),8,binsPlot);
  
  for (int i=0;i<etaLow.size();i++){
    histAsym->SetBinContent(i+2,asym[i]);
    histAsym->SetBinError(i+2,error[i]);
  }

  histAsym->SetTitle("");
  histAsym->SetXTitle("Lepton #eta");
  histAsym->SetYTitle("Lepton charge asymetry");  
  histAsym->GetYaxis()->CenterTitle();
  histAsym->GetXaxis()->CenterTitle();
  histAsym->SetMarkerStyle(20);
  if(doW) histAsym->GetYaxis()->SetRangeUser(0.02,0.34); 
  else histAsym->GetYaxis()->SetRangeUser(-0.1,0.1); 

  return histAsym;
}

// Function getAsymMC: computes the raw cut and count asymetry for a given MetCut and PtCut on MC
// ===============================================================================================

  
//////////////////////////////////////////////////////////////////////////////////////////
// TODO: write this part taking advantage of class CutAndCount and using example getAsym
///////////////////////////////////////////////////////////////////////////////////////////

TH1F* getAsymMC(float MetCut, float PtCut, bool doW){

  cout<<" This has to be coded ! "<< endl;

  TH1F *histAsym;  
  return histAsym;
}

// Function DrawCutAndCountAsym: Draws cut and count asymetries and compares to CMS paper results
// ===============================================================================================
void DrawCutAndCountAsym(bool doW){

  SetStyle();

  TH1F *Asym_25_25 =  getAsym(25,25, doW);
  TH1F *Asym_30_25 =  getAsym(30,25, doW);
  TH1F *Asym_25_30 =  getAsym(25,30, doW);
  TH1F *Asym_30_30 =  getAsym(30,30, doW);

  TH1F *Paper_30  =  getPaperPlot(30);
  TH1F *Paper_25  =  getPaperPlot(25);

  Asym_25_25->SetMarkerColor(kRed);
  Asym_25_25->SetLineColor(kRed);
  Asym_25_30->SetMarkerColor(kBlue);
  Asym_25_30->SetLineColor(kBlue);
  Paper_25->SetMarkerColor(kBlack);
  Paper_25->SetMarkerColor(kBlack);

  Asym_30_25->SetMarkerColor(kRed);
  Asym_30_25->SetLineColor(kRed);
  Asym_30_30->SetMarkerColor(kBlue);
  Asym_30_30->SetLineColor(kBlue);
  Paper_30->SetMarkerColor(kBlack);
  Paper_30->SetMarkerColor(kBlack);

  TLegend *leg25= new TLegend(0.2,0.7,0.5,0.85); 
  leg25->SetBorderSize(0);
  leg25->SetFillStyle(0);
  leg25->SetTextFont(42);
  leg25->SetHeader("pT > 25 GeV");
  if(doW){ 
    leg25->AddEntry(Asym_25_25,"M_{ET} > 25 GeV","lm"); 
    leg25->AddEntry(Asym_25_30,"M_{ET} > 30 GeV","lm"); 
    leg25->AddEntry(Paper_25,"CMS paper","lm"); 
  }
  TLegend *leg30= new TLegend(0.2,0.7,0.5,0.85); 
  leg30->SetBorderSize(0);
  leg30->SetFillStyle(0);
  leg30->SetTextFont(42);
  leg30->SetHeader("pT > 30 GeV");
  if(doW){ 
    leg30->AddEntry(Asym_30_25,"M_{ET} > 25 GeV","lm"); 
    leg30->AddEntry(Asym_30_30,"M_{ET} > 30 GeV","lm"); 
    leg30->AddEntry(Paper_30,"CMS paper","lm"); 
  }
  TCanvas *can25 = new TCanvas("can25","can25");
  Asym_25_25->Draw("E1");
  if(doW) Asym_25_30->Draw("E1same");
  if(doW) Paper_25->Draw("E1same");
  if(doW) leg25->Draw("same");
  if(doW)
    can25->Print("plots/WAsym_CutAndCount_Pt25.pdf","pdf");
  else
    can25->Print("plots/ZAsym_CutAndCount.pdf","pdf");

  TCanvas *can30 = new TCanvas("can30","can30");
  Asym_30_25->Draw("E1");
  if(doW) Asym_30_30->Draw("E1same");
  if(doW) Paper_30->Draw("E1same");
  if(doW)leg30->Draw("same");
  if(doW)
    can30->Print("plots/WAsym_CutAndCount_Pt30.pdf","pdf");
  else  
    can30->Print("plots/ZAsym_CutAndCount.pdf","pdf");
 
}



// Function DrawCutAndCountAsymMC: Draws cut and count asymetries on MC and compares to MC result no met cut
// ===========================================================================================================
void DrawCutAndCountAsymMC(){
  
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  // TODO: write this part taking advantage of example DrawCutAndCountAsym to draw W MC asymetries for different MET cuts and pt cuts
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  SetStyle();

  cout<<" This has to be coded ! "<< endl;

  bool doW=true;
  
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Exercice 2.1: Determine MET smearing to be applied to MC on a Z sample and determine PDFs for 3 species
///////////////////////////////////////////////////////////////////////////////////////////////////////////

// Function fitMetSmearing: (2.1.1) Determines the best smearing for a given ptCut
// ===============================================================================




void fitMetSmearing(float ptCut=25){

  SetStyle();
  // Met Histos for DY data:
  METSmearing *A = new METSmearing();
  TChain* ch = new TChain("events");
  ch->Add("files/data.root"); 
  A->SetSmearing(0);
  A->SetPtCut(ptCut);
  A->SetNMuons(2);
  ch->Process(A);
  TH1F *METData2Mu=(TH1F*) A->histograms.at(0)->Clone("METData2Mu");
  METData2Mu->SetLineWidth(3);
  METData2Mu->SetXTitle("M_{ET}");

  // Met Histos for DY MC:
  TChain* ch2 = new TChain("events");
  ch2->Add("files/dy.root"); 

  TH1F *METMC2Mu[30];
  TH1F *METMCSmeared2Mu[30];
  int nbins=30; //10
  float min=2.0; //5.0
  float max=8.0; // 6.0
  
  for(int i=0;i<nbins;i++){
    float smear=min+float(i)*(max-min)/nbins;
    METSmearing *B = new METSmearing();
    B->SetNMuons(2);
    B->SetPtCut(ptCut);
    B->SetSmearing(smear);
    ch2->Process(B);
    
    TString Name="METMCSmeared2Mu";
    Name+=i;
    METMCSmeared2Mu[i]=(TH1F*) B->histogramsSmeared.at(0)->Clone(Name.Data());
    METMCSmeared2Mu[i]->Scale(METData2Mu->Integral()/METMCSmeared2Mu[i]->Integral());
    METMCSmeared2Mu[i]->SetLineColor(i+2);    
    METMCSmeared2Mu[i]->SetXTitle("M_{ET}");
    
    METMC2Mu[i]=(TH1F*) B->histograms.at(0)->Clone(Name.Data());
    METMC2Mu[i]->Scale(METData2Mu->Integral()/METMC2Mu[i]->Integral());
    METMC2Mu[i]->SetLineColor(i+2);    
    METMC2Mu[i]->SetXTitle("M_{ET}");
  }

  

  TLegend *legNoSmear= new TLegend(0.5,0.7,0.8,0.85); 
  legNoSmear->SetBorderSize(0);
  legNoSmear->SetFillStyle(0);
  legNoSmear->SetTextFont(42);
  TString ls="pT >"; ls+=int(ptCut);ls+=" GeV";
  legNoSmear->SetHeader(ls.Data());
  legNoSmear->AddEntry(METData2Mu,"Z data","lm"); 
  legNoSmear->AddEntry( METMC2Mu[0],"Z MC No Smearing","f"); 

  TCanvas *canNoSmear = new TCanvas("canNoSmear","canNoSmear");
  METData2Mu->SetLineColor(kBlack);
  // METData2Mu->Draw();
  METMC2Mu[0]->SetFillColor(kRed);
  METMC2Mu[0]->Draw("H");
  METData2Mu->Draw("same");
  legNoSmear->Draw("same");
  if(ptCut==25) canNoSmear->Print("plots/noSmearFitSmearing25.pdf","pdf");
  else canNoSmear->Print("plots/noSmearFitSmearing30.pdf","pdf");

  TCanvas *can = new TCanvas("can","can");
  METData2Mu->SetLineColor(kBlack);
  METData2Mu->Draw();
  
  TH1D* chi2Hist = new TH1D("chi2Hist","chi2Hist",nbins,min-0.5*(max-min)/float(nbins),max-0.5*(max-min)/float(nbins));
  chi2Hist->SetXTitle("smearing");
  chi2Hist->SetYTitle("chi2");

  float chi2min=999;
  float bestSmear=0;
  int bestI=0;
  for(int i=0;i<nbins;i++){
    can->cd();
    METMCSmeared2Mu[i]->Draw("same");
    float smear=min+float(i)*(max-min)/nbins;
    cout<<"smear="<<smear<<" chi2test:"<<	METMCSmeared2Mu[i]->Chi2Test(METData2Mu,"NORMUUCHI2/NDF") << " prob:"<<	METMCSmeared2Mu[i]->Chi2Test(METData2Mu,"NORMUU") << " kolmogorov:"<< METMCSmeared2Mu[i]->KolmogorovTest(METData2Mu)<<endl; 
    float chi2=METMCSmeared2Mu[i]->Chi2Test(METData2Mu,"NORMUUCHI2/NDF");
    if(chi2<chi2min){
      chi2min=chi2;
      bestSmear=smear;
      bestI=i;
    }
    chi2Hist->SetBinContent(i+1,chi2);
    can->Update(); 
  }    

  cout<<" BEST SMEARING VALUE = "<<bestSmear<<" with chi2 = "<<chi2min<< endl;
  
  if(ptCut==25) can->Print("plots/fitSmearing25.pdf","pdf");
  else can->Print("plots/fitSmearing30.pdf","pdf");

  TLegend *legBest= new TLegend(0.5,0.7,0.8,0.85); 
  legBest->SetBorderSize(0);
  legBest->SetFillStyle(0);
  legBest->SetTextFont(42);
  TString ls2="pT >"; ls2+=int(ptCut);ls2+=" GeV";
  legBest->SetHeader(ls2.Data());
  legBest->AddEntry(METData2Mu,"Z data","lm"); 
  legBest->AddEntry( METMC2Mu[0],"Z MC Smeared","f"); 

  TCanvas *canBest = new TCanvas("canBest","canBest");
  METData2Mu->SetLineColor(kBlack);
  METData2Mu->Draw();
  METMCSmeared2Mu[bestI]->SetFillColor(kRed);
  METMCSmeared2Mu[bestI]->Draw("Hsame");
  METData2Mu->Draw("same");
  legBest->Draw("same");
  if(ptCut==25) canBest->Print("plots/bestFitSmearing25.pdf","pdf");
  else canBest->Print("plots/bestFitSmearing30.pdf","pdf");


  TCanvas*canChi2 = new TCanvas("canChi2","canChi2");
  chi2Hist->Draw("H");
  if(ptCut==25) canChi2->Print("plots/chi2FitSmearing25.pdf","pdf");
  else canChi2->Print("plots/chi2FitSmearing30.pdf","pdf");

}

// Function applyMetSmearing: (2.1.2) and (2.1.3) Applies best smearing and saves PDFs
// ===================================================================================

void applyMetSmearing(float ptCut){


  SetStyle();

  TFile *PDF;
  if(ptCut==25) PDF = new TFile("results/PDF.root","RECREATE");
  else PDF = new TFile("results/PDF30.root","RECREATE");
  
  // Met Histos for W data:
  METSmearing *A = new METSmearing();
  TChain* ch = new TChain("events");
  ch->Add("files/data.root"); 
  A->SetSmearing(0);
  A->SetPtCut(ptCut);
  A->SetNMuons(1);
  ch->Process(A);
  TH1F *METData1Mu=(TH1F*) A->histograms.at(0)->Clone("METData1Mu");
  METData1Mu->SetXTitle("M_{ET}");

  // best smearing:
  float smearingFromZ=5.7;
  if(ptCut==30) smearingFromZ=5.3;


  // Met Histos for W MC and smeared MET for W MC:
  METSmearing *B = new METSmearing();
  TChain* ch2 = new TChain("events");
  ch2->Add("files/wjets.root"); 
  B->SetSmearing(smearingFromZ);
  B->SetPtCut(ptCut);
  B->SetNMuons(1);
  ch2->Process(B);
  TH1F *METW1Mu=(TH1F*) B->histograms.at(0)->Clone("METW1Mu");
  METW1Mu->SetXTitle("M_{ET}");
  TH1F *METWSmeared1Mu=(TH1F*) B->histogramsSmeared.at(0)->Clone("PDF_W");
  METWSmeared1Mu->SetXTitle("M_{ET}");

  
  // Plot W:

  TLegend *leg= new TLegend(0.5,0.7,0.8,0.85); 
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  leg->SetTextFont(42);
  TString ls="pT >"; ls+=int(ptCut);ls+=" GeV";
  leg->SetHeader(ls.Data());
  leg->AddEntry(METData1Mu,"Data 1 #mu","lm"); 
  leg->AddEntry(METW1Mu ,"W MC","lm"); 
  leg->AddEntry(METWSmeared1Mu,"W MC Smeared","lm");

  TCanvas *can = new TCanvas("can","can");
  METData1Mu->SetLineColor(kBlack);
  METW1Mu->SetLineColor(kBlue);
  METW1Mu->Draw();
  METData1Mu->Draw("same");
  METWSmeared1Mu->SetLineColor(kRed);
  METWSmeared1Mu->Draw("same");
  leg->Draw("same");
  TString canname="plots/applySmearingToW";
  canname+=int(ptCut);
  canname+=".pdf";
  can->Print(canname.Data(),"pdf");
  

  // Met Histos for Z data with 1 muon and smeared histo:

  METSmearing *C = new METSmearing();
  TChain* ch3 = new TChain("events");
  ch3->Add("files/dy.root"); 
  C->SetSmearing(smearingFromZ);
  C->SetPtCut(ptCut);
  C->SetNMuons(1);
  ch3->Process(C);
  TH1F *METZ1Mu=(TH1F*) C->histograms.at(0)->Clone("METZ1Mu");
  TH1F *METZSmeared1Mu=(TH1F*) C->histogramsSmeared.at(0)->Clone("PDF_Z");
  METZ1Mu->SetXTitle("M_{ET}");
  METZSmeared1Mu->SetXTitle("M_{ET}");
  // Plot Z:

  TLegend *leg2= new TLegend(0.5,0.7,0.8,0.85); 
  leg2->SetBorderSize(0);
  leg2->SetFillStyle(0);
  leg2->SetTextFont(42);
  leg2->SetHeader(ls.Data());
  //leg2->AddEntry(METData1Mu,"Data 1 #mu","lm"); 
  leg2->AddEntry(METZ1Mu ,"Z MC 1 #mu","lm"); 
  leg2->AddEntry(METWSmeared1Mu,"Z MC 1 #mu Smeared","lm");

  TCanvas *can2 = new TCanvas("can2","can2");
  METData1Mu->SetLineColor(kBlack);
  METZ1Mu->SetLineColor(kBlue);
  METZSmeared1Mu->SetLineColor(kRed);
  METZ1Mu->Draw();
  METZSmeared1Mu->Draw("same");
  leg2->Draw("same");
  TString can2name="plots/applySmearingToZ";
  can2name+=int(ptCut);
  can2name+=".pdf";
  can2->Print(can2name.Data(),"pdf");


  // Met Histos for data with 0 muon:

  METSmearing *D = new METSmearing();
  TChain* ch4 = new TChain("events");
  ch4->Add("files/data.root"); 
  D->SetSmearing(0);
  D->SetPtCut(ptCut);
  D->SetNMuons(0);
  ch4->Process(D);
  TH1F *METData0Mu=(TH1F*) D->histograms.at(0)->Clone("METData0Mu");  
  METData0Mu->SetXTitle("M_{ET}");

  // Met Histos for W MC with 0 muon:
  METSmearing *E = new METSmearing();
  TChain* ch5 = new TChain("events");
  ch5->Add("files/wjets.root"); 
  E->SetSmearing(0);
  E->SetPtCut(ptCut);
  E->SetNMuons(0);
  ch5->Process(E);
  TH1F *METW0Mu=(TH1F*) E->histograms.at(0)->Clone("METW0Mu");
  METW0Mu->SetXTitle("M_{ET}");

  // Compute QCD PDF data(0mu)-W MC(0mu)
  TH1F *PDF_QCD=new TH1F("PDF_QCD","PDF_QCD",METW0Mu->GetNbinsX(),METW0Mu->GetBinLowEdge(1),METW0Mu->GetBinLowEdge(METW0Mu->GetNbinsX())+METW0Mu->GetBinWidth(METW0Mu->GetNbinsX()));
  for(int j=0;j<METW0Mu->GetNbinsX();j++){
    float bc=METData0Mu->GetBinContent(j+1)-METW0Mu->GetBinContent(j+1); float be=METData0Mu->GetBinError(j+1);
    if(bc<0) bc=0;
    PDF_QCD->SetBinContent(j+1,bc);
    PDF_QCD->SetBinError(j+1,be);
  }
  PDF_QCD->SetXTitle("M_{ET}");

  // Draw QCD:
  TLegend *leg3= new TLegend(0.5,0.7,0.8,0.85); 
  leg3->SetBorderSize(0);
  leg3->SetFillStyle(0);
  leg3->SetTextFont(42);
  leg3->SetHeader(ls.Data());
  leg3->AddEntry(METData0Mu,"Data 0 #mu","lm"); 
  leg3->AddEntry(METW0Mu ,"W MC 0 #mu","lm"); 
  leg3->AddEntry(PDF_QCD,"QCD PDF","lm");

  TCanvas *can3 = new TCanvas("can3","can3");
  METData0Mu->SetLineColor(kBlack);
  METW0Mu->SetLineColor(kBlue);
  PDF_QCD->SetLineColor(kRed);
  METData0Mu->Draw();
  METW0Mu->Draw("same");
  PDF_QCD->Draw("same"); 
  leg3->Draw("same");
  TString can3name="plots/QCD";
  can3name+=int(ptCut);
  can3name+=".pdf";
  can3->Print(can3name.Data(),"pdf");

  METZSmeared1Mu->Write();
  METWSmeared1Mu->Write();
  PDF_QCD->Write();
  PDF->Close();

}


////////////////////////////////////////
// Exercice 2.2: Maximum Likelihood Fit
////////////////////////////////////////

// Function getZfrac: (2.2.1) Determines the ratio between the total Z and W yields versus eta 
// ============================================================================================



void getZfrac(float PtCut){

  SetStyle();
  // Compute W numbers for a given pt cut 
  CutAndCount *A = new CutAndCount();
  TChain* ch = new TChain("events");
  ch->Add("files/wjets.root"); 
  A->SetPtCut(PtCut);
  A->SetMetCut(0);
  A->SetDoW(true);
  ch->Process(A);

  vector<float> nWplus=A->nWplus;
  vector<float> nWminus=A->nWminus;

  // Compute Z numbers for a given pt cut
  CutAndCount *B = new CutAndCount();
  TChain* ch2 = new TChain("events");
  ch2->Add("files/dy.root"); 
  B->SetPtCut(PtCut);
  B->SetMetCut(0);
  B->SetDoW(true);
  ch2->Process(B);
  
  vector<float> nZplus=B->nWplus;
  vector<float> nZminus=B->nWminus;
  
  vector<float> etaLow=A->EtaLow;
  vector<float> etaHigh=A->EtaHigh;
 
  int nbins=6;
  float binsPlot[9];
  binsPlot[0]=-0.1;
  binsPlot[8]=2.7;

  for (int i=0;i<etaLow.size();i++){

    binsPlot[i+1]=etaLow[i];
  }
  
  binsPlot[nbins+1]=etaHigh[nbins-1];

  
  TString nameh="Zfrac_";
  nameh+=PtCut;

  TH1F *histZfrac= new TH1F(nameh.Data(),nameh.Data(),8,binsPlot);

  for (int j=0;j<6;j++){
    histZfrac->SetBinContent(j+2,(nZplus[j]+nZminus[j])/(nWplus[j]+nWminus[j]));
    histZfrac->SetBinError(j+2,((nZplus[j]+nZminus[j])/(nWplus[j]+nWminus[j]))*sqrt((1./(nWplus[j]+nWminus[j]))+(1./(nZplus[j]+nZminus[j]))));
  }

  TCanvas *can = new TCanvas("can","can");

  histZfrac->SetXTitle("Lepton #eta");
  histZfrac->SetYTitle("n_{Z}/n_{W}");  
  histZfrac->SetMarkerStyle(21);
  histZfrac->GetYaxis()->CenterTitle();
  histZfrac->GetXaxis()->CenterTitle();
  histZfrac->Draw("E1");
  TString nameh2="plots/Zfrac_";
  nameh2+=PtCut;
  nameh2+=".pdf";
  can->Print(nameh2.Data(),"pdf");

  // Compute ratio:
  for(int j=0;j<nZplus.size();j++){  
    cout<<" Z frac"<<  A->EtaLow[j]<<" < eta < "<<A->EtaHigh[j]<< " = "<<(nZplus[j]+nZminus[j])/(nWplus[j]+nWminus[j])<<" +/- "<< ((nZplus[j]+nZminus[j])/(nWplus[j]+nWminus[j]))*sqrt((1./(nWplus[j]+nWminus[j]))+(1./(nZplus[j]+nZminus[j])))<< endl;
  }
}


void DrawPDF(float ptCut){

  SetStyle();
  TString ls="pT >"; ls+=int(ptCut);ls+=" GeV";

  TFile *file;
  if(ptCut==25) file= new TFile("results/PDF.root");
  else file= new TFile("results/PDF30.root");
  TH1F *PDF_Z=(TH1F*)file->Get("PDF_Z");
  TH1F *PDF_W=(TH1F*)file->Get("PDF_W");
  TH1F *PDF_QCD=(TH1F*)file->Get("PDF_QCD");
  PDF_Z->Scale(1./PDF_Z->Integral());
  PDF_W->Scale(1./PDF_W->Integral());
  PDF_QCD->Scale(1./PDF_QCD->Integral());
  PDF_W->SetLineColor(kRed);
  PDF_W->SetLineWidth(2);
  PDF_Z->SetLineColor(kBlue);
  PDF_Z->SetLineWidth(2);
  PDF_QCD->SetLineColor(kBlack);
  PDF_QCD->SetLineWidth(2);


  TLegend *leg= new TLegend(0.5,0.7,0.8,0.85); 
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  leg->SetTextFont(42);
  leg->SetHeader(ls.Data());
  leg->AddEntry(PDF_W,"PDF W","lm");
  leg->AddEntry(PDF_QCD,"PDF QCD","lm"); 
  leg->AddEntry(PDF_Z ,"PDF Z","lm"); 

  TCanvas *can = new TCanvas("can","can");
  PDF_QCD->Draw();
  PDF_W->Draw("same");
  PDF_Z->Draw("same");leg->Draw("same");
  if(ptCut==25) can->Print("plots/PDFs.pdf","pdf");
  else can->Print("plots/PDFs30.pdf","pdf");
}





// Function makeTree: to make a simpler data tree for the fit (tree already done)
// ==============================================================================

void makeTree(){

  MakeTree *A = new MakeTree();
  TChain* ch = new TChain("events");
  ch->Add("files/data.root");
  ch->Process(A);
}

// NB: the rest of exercice 2 will be done in a separate macro.


int main() {

  // To draw CMS paper plot:
  //=========================
  //DrawPaperPlot();

  // Exercice 1.1:
  //===============
  //BasicStackPlots();

  // Exercice 1.2:
  //===============
  //DrawCutAndCountAsym(true);
  //DrawCutAndCountAsymMC();

  // Exercice 2.1:
  //===============

  fitMetSmearing(25);  
  //fitMetSmearing(30); 

  applyMetSmearing(25);
  //applyMetSmearing(30);

  DrawPDF(25);
  //DrawPDF(30);
  
  // Exercice 2.2:
  //===============
  //getZfrac(25);
  //getZfrac(30);  

}