//////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #202
// 
// Setting up an extended maximum likelihood fit
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooDataHist.h"
#include "RooHistPdf.h"
#include "RooFitResult.h"
#include "RooExtendPdf.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TStyle.h"
#include "TH1F.h"
#include "TLegend.h"
#include <TROOT.h>
#include "TAxis.h"
#include <vector>
#include <TGraphAsymmErrors.h>
#include "RooPlot.h"
using namespace RooFit ;
using namespace std;


void SetStyle(){

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

  // gStyle->SetTitleOffset(0.7,"x");
  //gStyle->SetTitleOffset(0.7,"y");

  gStyle->SetPadBottomMargin(0.15);
  gStyle->SetPadTopMargin(0.05);
  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadRightMargin(0.05);
  
  gStyle->SetTitleColor(1);
  gStyle->SetTitleFillColor(0);
  gStyle->SetTitleFontSize(0.05);
  gStyle->SetTitleBorderSize(0);
  
}

vector<float> simFitMET( float ptCut, int etaBin, bool zFracFree)
{
  gROOT->Reset();
  SetStyle();

  // bool zFracFree=false;
   vector<float> etabounds;
   etabounds.push_back(0);
   etabounds.push_back(0.4);
   etabounds.push_back(0.8);
   etabounds.push_back(1.2);
   etabounds.push_back(1.5);
   etabounds.push_back(1.8);
   etabounds.push_back(2.1);

   vector<float> zfrac;
   vector<float> zfracLow;
   vector<float> zfracHigh;

   if(ptCut==25.){

     zfrac.push_back(0.046);
     zfrac.push_back(0.051);
     zfrac.push_back(0.066);
     zfrac.push_back(0.076);
     zfrac.push_back(0.085);
     zfrac.push_back(0.096);

     zfracLow.push_back(0.0424114);
     zfracHigh.push_back(0.0514919);
     zfracLow.push_back(0.0474304);
     zfracHigh.push_back(0.0557186);
     zfracLow.push_back(0.0589267);
     zfracHigh.push_back(0.0767653);
     zfracLow.push_back(0.0677544);
     zfracHigh.push_back(0.0876032);
     zfracLow.push_back(0.0703208); 
     zfracHigh.push_back(0.10801); 
     zfracLow.push_back(0.0809825);
     zfracHigh.push_back(0.1228);

   }else{

     zfrac.push_back(0.039);
     zfrac.push_back(0.042);
     zfrac.push_back(0.059);
     zfrac.push_back(0.070);
     zfrac.push_back(0.082);
     zfrac.push_back(0.097);

     zfracLow.push_back(0.0353328);
     zfracHigh.push_back(0.0433999);
     zfracLow.push_back(0.0395586);
     zfracHigh.push_back(0.0446883);
     zfracLow.push_back(0.052955);
     zfracHigh.push_back(0.0681026);
     zfracLow.push_back(0.0638405);
     zfracHigh.push_back(0.078731);
     zfracLow.push_back(0.0687565);
     zfracHigh.push_back(0.10319);
     zfracLow.push_back(0.0830126 );
     zfracHigh.push_back(0.120892);

   }

   // S e t u p   c o m p o n e n t   p d f s 
  // ---------------------------------------

  // define dataset to fit

  // Define observables y,z
  RooRealVar met("MET","M_{ET}",0,120.) ;
  RooRealVar charge("charge","charge",-1,1);
  RooRealVar pt("pt","pt",0,300.0);
  RooRealVar eta("eta","eta",0,3.0);

  // Import only observables (y,z)
  TFile *fileData = new TFile("files/dataForFit.root");
  TTree *treeData = (TTree*) fileData->Get("data");

  RooDataSet dataAll("dataAll","dataAll",RooArgSet(met,charge,pt,eta),Import(*treeData)) ;  
  TString cut=" eta>= ";
  cut+=etabounds[etaBin];
  cut+=" && eta< ";
  cut+=etabounds[etaBin+1];
  cut+=" && pt>";
  cut+=ptCut;  
  cout<<" CUT="<< cut<< endl;
  TString cutPlus="charge == 1 &&";
  cutPlus+=cut;
  TString cutMinus="charge == -1 &&";
  cutMinus+=cut;

  RooDataSet *dataPlus= (RooDataSet*) dataAll.reduce(cutPlus.Data());
  RooDataSet *dataMinus= (RooDataSet*) dataAll.reduce(cutMinus.Data());
  dataPlus->Print() ;
  dataMinus->Print() ;

  // C r e a t e   i n d e x   c a t e g o r y   a n d   j o i n   s a m p l e s 
  // ---------------------------------------------------------------------------

  // Define category to distinguish physics and control samples events
  RooCategory WCharge("WCharge","WCharge") ;
  WCharge.defineType("plus") ;
  WCharge.defineType("minus") ;

  // Construct combined dataset in (x,sample)
  RooDataSet combData("combData","combined data",met,Index(WCharge),Import("plus",*dataPlus),Import("minus",*dataMinus)) ;

  // open pdfs file
  TFile *filePdf;
  if(ptCut==25) filePdf = new TFile("results/PDF.root");
  else filePdf = new TFile("results/PDF30.root");
  TH1F *PDF_Z = (TH1F*) filePdf->Get("PDF_Z");
  TH1F *PDF_W = (TH1F*) filePdf->Get("PDF_W");
  TH1F *PDF_QCD = (TH1F*) filePdf->Get("PDF_QCD");
  
  RooDataHist *histZ = new RooDataHist("histZ","histZ", RooArgList(met), PDF_Z);  
  RooDataHist *histW   = new RooDataHist("histW",  "histW",   RooArgList(met), PDF_W);
  RooDataHist *histQCD = new RooDataHist("histQCD","histQCD", RooArgList(met), PDF_QCD); 

  // Represent data in dh as pdf in x
  RooHistPdf pdfZ("pdfZ","pdfZ",met,*histZ,0) ;
  RooHistPdf pdfW("pdfW","pdfW",met,*histW,0) ;
  RooHistPdf pdfQCD("pdfQCD","pdfQCD",met,*histQCD,0) ;

  // Sum the signal components into a composite signal p.d.f.

  RooRealVar nW("nW","Total Signal Yield",100,0.,211000.) ;
  RooRealVar Asym("Asym","W charge asymetry",0.2,0.,1.) ;  
  RooRealVar Zfrac("Zfrac","fraction of Z wrt W", zfrac[etaBin],zfracLow[etaBin],zfracHigh[etaBin]);
  if(!zFracFree) Zfrac.setConstant(kTRUE);


  RooFormulaVar nZ("nZ","Total Z Yield","@0*@1",RooArgList(Zfrac,nW));
  RooRealVar nQCD("nQCD","QCD Yield",100,0.,2000000.);
  RooRealVar half("half","half",0.5);

  RooFormulaVar nWPlus("nWPlus","Plus Signal Yield","@0*0.5*(1+@1)", RooArgList(nW,Asym));
  RooFormulaVar nWMinus("nWMinus","Minus Signal Yield","@0*0.5*(1-@1)", RooArgList(nW,Asym));
  RooFormulaVar nZPoM("nZPoM","Half Z Yield","@0*@1",RooArgList(nZ,half));
  RooFormulaVar nQCDPoM("nQCDPoM","HalfQCDYield","@0*@1",RooArgList(nQCD,half));

  /////////////////////
  // M E T H O D   1 //
  /////////////////////

  
  // C o n s t r u c t   e x t e n d e d   c o m p o s i t e   m o d e l 
  // -------------------------------------------------------------------
  
  RooAddPdf  modelPlus("modelPlus","modelPlus",RooArgList(pdfW,pdfZ,pdfQCD),RooArgList(nWPlus,nZPoM,nQCDPoM)) ;
  RooAddPdf  modelMinus("modelMinus","modelMinus",RooArgList(pdfW,pdfZ,pdfQCD),RooArgList(nWMinus,nZPoM,nQCDPoM)) ;
  
  // Construct a simultaneous pdf using category sample as index
  RooSimultaneous simPdf("simPdf","simultaneous pdf",WCharge) ;

  // Associate model with the physics state and model_ctl with the control state
  simPdf.addPdf(modelPlus,"plus") ;
  simPdf.addPdf(modelMinus,"minus") ;


  // S a m p l e ,   f i t   a n d   p l o t   e x t e n d e d   m o d e l 
  // ---------------------------------------------------------------------

  // Fit model to data, extended ML term automatically included
  RooFitResult*fitRes=simPdf.fitTo(combData,Save(),Minos(1)) ; 
  
  // Plot data and PDF overlaid, use expected number of events for p.d.f projection normalization
  // rather than observed number of events (==data->numEntries())
  
  RooPlot* xframe = met.frame(Title("W plus")) ;
  dataPlus->plotOn(xframe, Name("data")) ;
  modelPlus.plotOn(xframe,Name("model"),Normalization(1.0,RooAbsReal::RelativeExpected)) ;

  // Overlay the background component of model with a dashed line
  modelPlus.plotOn(xframe,Components(pdfZ),Name("Z"),LineStyle(kDashed),LineColor(kGreen),Normalization(1.0,RooAbsReal::RelativeExpected)) ;

  // Overlay the background+sig2 components of model with a dotted line
  modelPlus.plotOn(xframe,Components(RooArgSet(pdfZ,pdfQCD)),Name("ZQCD"),LineStyle(kDotted),LineColor(kRed),Normalization(1.0,RooAbsReal::RelativeExpected)) ;


  //legPlus->AddEntry(histAsym,"Result","lm"); 
  //legPlus->AddEntry(Paper_25,"CMS paper","lm"); 
  //legPlus->AddEntry(theo25[0],"Resbos","fl"); 
  //legPlus->AddEntry(theo25[1],"MCFM","fl"); 
  
  RooPlot* xframe2 = met.frame(Title("W minus")) ;
  dataMinus->plotOn(xframe2) ;
  modelMinus.plotOn(xframe2,Name("model"),Normalization(1.0,RooAbsReal::RelativeExpected)) ;

  // Overlay the background component of model with a dashed line
  modelMinus.plotOn(xframe2,Components(pdfZ),Name("Z"),LineStyle(kDashed),LineColor(kGreen),Normalization(1.0,RooAbsReal::RelativeExpected)) ;

  // Overlay the background+sig2 components of model with a dotted line
  modelMinus.plotOn(xframe2,Components(RooArgSet(pdfZ,pdfQCD)),Name("ZQCD"),LineStyle(kDotted),LineColor(kRed),Normalization(1.0,RooAbsReal::RelativeExpected)) ;

  // Print structure of composite p.d.f.
  //model.Print("t") ;

  



  // Draw the frame on the canvas
  TCanvas *canPlus = new TCanvas("canPlus","canPlus",600,600) ;
  gPad->SetLeftMargin(0.15) ; gPad->SetBottomMargin(0.15) ; xframe->GetYaxis()->SetTitleOffset(1.) ; xframe->GetXaxis()->SetTitleOffset(1.0) ;  xframe->Draw() ;
  
  TLegend *legPlus= new TLegend(0.55,0.5,0.9,0.9); 
  legPlus->SetBorderSize(0);
  legPlus->SetFillStyle(0);
  legPlus->SetTextFont(42);
  TString ptPlus="W+, pT >";
  ptPlus+=int(ptCut);
  ptPlus+= " GeV , #eta bin ";
  ptPlus+=etaBin;
  legPlus->SetHeader(ptPlus.Data());   
  legPlus->AddEntry("data","Data", "PLM");
  legPlus->AddEntry(xframe->findObject("Z"),"Z background", "L");
  legPlus->AddEntry(xframe->findObject("ZQCD"),"Z + QCD backgrounds", "L");
  legPlus->AddEntry(xframe->findObject("model"),"Signal + backgrounds", "L");
  legPlus->Draw();
  
  TCanvas *canMinus = new TCanvas("canMinus","canMinus",600,600) ;
  gPad->SetLeftMargin(0.15) ; gPad->SetBottomMargin(0.15) ; xframe2->GetYaxis()->SetTitleOffset(1.) ; xframe->GetXaxis()->SetTitleOffset(1.0) ; xframe2->Draw() ; 
  TLegend *legMinus= new TLegend(0.55,0.5,0.9,0.9); 
  legMinus->SetBorderSize(0);
  legMinus->SetFillStyle(0);
  legMinus->SetTextFont(42);
  TString ptMinus="W -, pT >";
  ptMinus+=int(ptCut);
  ptMinus+= " GeV , #eta bin ";
  ptMinus+=etaBin;
  legMinus->SetHeader(ptMinus.Data());   
  legMinus->AddEntry("data","Data", "PLM");
  legMinus->AddEntry(xframe2->findObject("Z"),"Z background", "L");
  legMinus->AddEntry(xframe2->findObject("ZQCD"),"Z + QCD backgrounds", "L");
  legMinus->AddEntry(xframe2->findObject("model"),"Signal + backgrounds", "L");
  legMinus->Draw();

  
  TString canNamePlus="plots/simFitMET_Plus_Bin";
  canNamePlus+=etaBin;
  canNamePlus+="Pt";
  canNamePlus+=int(ptCut);
  if(zFracFree) canNamePlus+="_zFracFree";
  canNamePlus+=".pdf";
  TString canNameMinus="plots/simFitMET_Minus_Bin";
  canNameMinus+=etaBin;
  canNameMinus+="Pt";
  canNameMinus+=int(ptCut);
  if(zFracFree) canNameMinus+="_zFracFree";
  canNameMinus+=".pdf";

  canPlus->Print(canNamePlus.Data());
  canMinus->Print(canNameMinus.Data());
  
   RooRealVar* AW = (RooRealVar*) fitRes->floatParsFinal().find("Asym");
   float asym=AW->getVal();
   float easym=AW->getError();
   vector<float> res;
   res.push_back(asym);
   res.push_back(easym);
   
   
   return res;
   
}


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;

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


void doFits(float ptCut, bool zFracFree){

  SetStyle();

   vector<float> etabounds;
   etabounds.push_back(0);
   etabounds.push_back(0.4);
   etabounds.push_back(0.8);
   etabounds.push_back(1.2);
   etabounds.push_back(1.5);
   etabounds.push_back(1.8);
   etabounds.push_back(2.1);

   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 values[6];
   float errors[6];
   
   for(int j=0;j<6;j++){
     vector<float> res=simFitMET(ptCut,j,zFracFree);
     values[j]=res[0];
     errors[j]=res[1];
   }
   
   TH1F *histAsym = new TH1F("Asym","Asym",8,binsPlot);
   
   for (int i=0;i<6;i++){
     histAsym->SetBinContent(i+2,values[i]);
     histAsym->SetBinError(i+2,errors[i]);
   }
   TH1F *Paper_25  =  getPaperPlot(ptCut);
   histAsym->SetMarkerColor(kRed);
   histAsym->SetLineColor(kRed);
   Paper_25->SetMarkerColor(kBlack);
   Paper_25->SetMarkerColor(kBlack);
   histAsym->SetTitle("");
   Paper_25->SetTitle("");
   
  vector<TGraphAsymmErrors*> theo25=getPaperPlotTheo(ptCut);

   histAsym->SetXTitle("Lepton #eta");
   histAsym->SetYTitle("Lepton charge asymetry");  

  histAsym->SetMarkerStyle(20);
  histAsym->GetYaxis()->SetRangeUser(0.02,0.34);
  histAsym->GetYaxis()->CenterTitle();
  histAsym->GetXaxis()->CenterTitle();

   TLegend *leg25= new TLegend(0.2,0.6,0.6,0.85); 
   leg25->SetBorderSize(0);
   leg25->SetFillStyle(0);
   leg25->SetTextFont(42);
   TString pt=" pT >";
   pt+=int(ptCut);
   pt+= " GeV";
   
   leg25->SetHeader(pt.Data());   
   leg25->AddEntry(histAsym,"Result","lm"); 
   leg25->AddEntry(Paper_25,"CMS paper","lm"); 
   leg25->AddEntry(theo25[0],"Resbos","fl"); 
   leg25->AddEntry(theo25[1],"MCFM","fl"); 
   
   TCanvas *can25 = new TCanvas("can25","can25");
   gPad->SetLeftMargin(0.15) ; gPad->SetBottomMargin(0.15) ; 
   histAsym->GetYaxis()->SetTitleOffset(1.1) ; 
   histAsym->GetXaxis()->SetTitleOffset(1.) ; 

   histAsym->Draw("E1");
   Paper_25->Draw("E1same");
   theo25[1]->Draw("3same");
   theo25[1]->Draw("Lsame");
   Paper_25->Draw("E1Same");
   histAsym->Draw("E1Same");
   leg25->Draw("same"); 
   TString cn="plots/WAsym_simFitMET_Pt";
   cn+=int(ptCut);
   if(zFracFree) cn+="_zFracFree";
   cn+=".pdf";
   can25->Print(cn.Data(),"pdf");
   
}

void test(){

  vector<float> res=simFitMET(25,0,false);

}