#include "TFile.h"
#include "TGraphAsymmErrors.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TAxis.h"
#include <iostream>
#include <sstream>
#include <iomanip>
#include "AtlasUtils.h"
#include "AtlasStyle.C"
#include "AtlasLabels.C"

#include "ComputeTotalCrossSection.C"
#include "readVectorsFromFile.C"

#define saveAsRoot      false
#define saveAsEps       true
#define showMCNLOnoMPI  false
#define m_debug         0
#define writeDataResults true

void makeDistribPlot(std::string hname="j_bjet_Pt_binned", bool setLogY=false, bool useOldaMCNLO=true)
{

  bool oldaMCNLO = useOldaMCNLO;
  double mpiGlobCorr = 1.;
  if (hname=="j_bjet_ZY_binned") {
      oldaMCNLO=false;
      mpiGlobCorr = 1.0286232;
  }
  
  gROOT->LoadMacro("AtlasUtils.C");
  SetAtlasStyle(); 

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


  Double_t Lumi = 4.6; // fb 

  if (hname=="XXX"){
    std::cout<<"No histogram selected - possible choices are: \n";
    std::cout<<"\"j_bjet_Pt_binned\""<<std::endl;
    std::cout<<"\"j_bjet_Y_binned\""<<std::endl;
    std::cout<<"\"j_bjet_ZPt_binned\""<<std::endl;
    std::cout<<"\"j_bjet_ZY_binned\""<<std::endl;
    std::cout<<"\"j_bjet_ZYf_binned\""<<std::endl;
    std::cout<<"\"j_bjet_ZdR620_binned\""<<std::endl;
    std::cout<<"\"j_bjet_ZdPhiF20_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_ZPt_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_dR_binned\""<<std::endl; 
    std::cout<<"\"j_2bjet_Mbb_binned\""<<std::endl;

    return;
  }

     

     TFile* fileSHPY = new TFile(  "../SherpaPythia/SherpaAndPythiaResultsNP.root");
     TFile* fileaMCNLOMPI; 
     if (oldaMCNLO) fileaMCNLOMPI = new TFile("../aMCatNLO/merged.root");
     else fileaMCNLOMPI = new TFile("../aMCatNLO/Zb.root");
     TFile* fileaMCNLO   = new TFile("../aMCatNLO/Zbb.root");
     TFile* fileMCFM;
     TFile* fileMCFMmstw;
     TFile* fileMCFMnnpdf;
     if (hname.substr(0,3)=="j_b") 
     {
       fileMCFM = new TFile("../MCFM/Zb/ZbMCFMgraphCT10.root");
       fileMCFMmstw = new TFile("../MCFMmstw/Zb/ZbMCFMgraphMSTW2008.root");
       fileMCFMnnpdf = new TFile("../MCFMnnpdf/Zb/ZbMCFMgraphNNPDF23.root");
     }
     else if (hname.substr(0,3)=="j_2")
     {
       fileMCFM = new TFile("../MCFM/Zbb/ZbbMCFMgraphCT10.root");       
       fileMCFMmstw = new TFile("../MCFMmstw/Zbb/ZbbMCFMgraphMSTW2008.root");
       fileMCFMnnpdf = new TFile("../MCFMnnpdf/Zbb/ZbbMCFMgraphNNPDF23.root");
     }



     TFile* foutRoot;
     if (saveAsRoot) 
       {
	 std::string fn = "MCFM_NP_FSR_"+hname+".root";
	 foutRoot = new TFile(fn.c_str(),"RECREATE");
       }


     // sherpa and pythia NP 
     TGraphAsymmErrors *grNPSherpa; 
     TGraphAsymmErrors *grNPPythia; 
     TGraphAsymmErrors *grNPstat;  
     TGraphAsymmErrors *grNPstatAndSyst;  

     // sherpa and pythia without NP 
     TGraphAsymmErrors *grSherpa; 
     TGraphAsymmErrors *grPythia; 
     TGraphAsymmErrors *graMCNLO; 
     TGraphAsymmErrors *grSherpafb; // in fb 
     TGraphAsymmErrors *grPythiafb; // in fb 
     TGraphAsymmErrors *graMCNLOfb; 
     // sherpa and pythia with NP 
     TGraphAsymmErrors *grSherpaNP; 
     TGraphAsymmErrors *grPythiaNP; 
     TGraphAsymmErrors *graMCNLOMPI; 
     TGraphAsymmErrors *grSherpaNPfb; // in fb 
     TGraphAsymmErrors *grPythiaNPfb; // in fb 
     TGraphAsymmErrors *graMCNLOMPIfb; // in fb 
     std::string hnameSherpa = "XS_sherpa_NPnone_"+hname;
     std::string hnamePythia = "XS_pythia_NPnone_"+hname;
     std::string hnameSherpaNP = "XS_sherpa_NPall_"+hname;
     std::string hnamePythiaNP = "XS_pythia_NPall_"+hname;
     std::string hnameNPSherpa = "NP_sherpa_"+hname;
     std::string hnameNPPythia = "NP_pythia_"+hname;

     TGraphAsymmErrors *grMCFM; 
     TGraphAsymmErrors *grMCFMas;   // with a_s error 
     TGraphAsymmErrors *grMCFMpdf;  // with pdf error 
     TGraphAsymmErrors *grMCFMsca;  // with scale error 
     TGraphAsymmErrors *grMCFMsta;  // with statistical error 
     // MCFM corrected for NP 
     TGraphAsymmErrors *grMCFMNPsta;  // with statistical error 
     TGraphAsymmErrors *grMCFMNPsta_as;   // AND with a_s error     
     TGraphAsymmErrors *grMCFMNPsta_as_pdf;  // AND with pdf error 
     TGraphAsymmErrors *grMCFMNPsta_as_pdf_sca;  // AND with scale error 
     TGraphAsymmErrors *grMCFMNPsta_as_pdf_sca_np;  // AND with np error 
     // MCFM with mstw2008 
     TGraphAsymmErrors *grMCFM8; 
     TGraphAsymmErrors *grMCFM8as;   // with a_s error 
     TGraphAsymmErrors *grMCFM8pdf;  // with pdf error 
     TGraphAsymmErrors *grMCFM8sca;  // with scale error 
     TGraphAsymmErrors *grMCFM8sta;  // with statistical error 
     // MCFM with mstw2008 corrected for NP 
     TGraphAsymmErrors *grMCFM8NPsta;  // with statistical error 
     TGraphAsymmErrors *grMCFM8NPsta_as;   // AND with a_s error     
     TGraphAsymmErrors *grMCFM8NPsta_as_pdf;  // AND with pdf error 
     TGraphAsymmErrors *grMCFM8NPsta_as_pdf_sca;  // AND with scale error 
     TGraphAsymmErrors *grMCFM8NPsta_as_pdf_sca_np;  // AND with np error 
     // MCFM with nnpdf23 
     TGraphAsymmErrors *grMCFM23; 
     TGraphAsymmErrors *grMCFM23as;   // with a_s error 
     TGraphAsymmErrors *grMCFM23pdf;  // with pdf error 
     TGraphAsymmErrors *grMCFM23sca;  // with scale error 
     TGraphAsymmErrors *grMCFM23sta;  // with statistical error 
     // MCFM with nnpdf23 corrected for NP 
     TGraphAsymmErrors *grMCFM23NPsta;  // with statistical error 
     TGraphAsymmErrors *grMCFM23NPsta_as;   // AND with a_s error     
     TGraphAsymmErrors *grMCFM23NPsta_as_pdf;  // AND with pdf error 
     TGraphAsymmErrors *grMCFM23NPsta_as_pdf_sca;  // AND with scale error 
     TGraphAsymmErrors *grMCFM23NPsta_as_pdf_sca_np;  // AND with np error 

     Double_t MCFMNP_Y[100];
     Double_t MCFMNP_EYlstat[100];
     Double_t MCFMNP_EYlstat_as[100];
     Double_t MCFMNP_EYlstat_as_pdf[100];
     Double_t MCFMNP_EYlstat_as_pdf_sca[100];
     Double_t MCFMNP_EYlstat_as_pdf_sca_np[100];
     Double_t MCFMNP_EYhstat[100];
     Double_t MCFMNP_EYhstat_as[100];
     Double_t MCFMNP_EYhstat_as_pdf[100];
     Double_t MCFMNP_EYhstat_as_pdf_sca[100];
     Double_t MCFMNP_EYhstat_as_pdf_sca_np[100];
     Double_t MCFM8NP_Y[100];
     Double_t MCFM8_EYlstat[100]; //// to be filled by scaling MCFM mstw2008
     Double_t MCFM8NP_EYlstat[100];
     Double_t MCFM8NP_EYlstat_as[100];
     Double_t MCFM8NP_EYlstat_as_pdf[100];
     Double_t MCFM8NP_EYlstat_as_pdf_sca[100];
     Double_t MCFM8NP_EYlstat_as_pdf_sca_np[100];
     Double_t MCFM8_EYhstat[100]; //// to be filled by scaling MCFM mstw2008
     Double_t MCFM8NP_EYhstat[100];
     Double_t MCFM8NP_EYhstat_as[100];
     Double_t MCFM8NP_EYhstat_as_pdf[100];
     Double_t MCFM8NP_EYhstat_as_pdf_sca[100];
     Double_t MCFM8NP_EYhstat_as_pdf_sca_np[100];
     Double_t MCFM23NP_Y[100];
     Double_t MCFM23_EYlstat[100]; //// to be filled by scaling MCFM nnpdf23
     Double_t MCFM23NP_EYlstat[100];
     Double_t MCFM23NP_EYlstat_as[100];
     Double_t MCFM23NP_EYlstat_as_pdf[100];
     Double_t MCFM23NP_EYlstat_as_pdf_sca[100];
     Double_t MCFM23NP_EYlstat_as_pdf_sca_np[100];
     Double_t MCFM23_EYhstat[100]; //// to be filled by scaling MCFM nnpdf23
     Double_t MCFM23NP_EYhstat[100];
     Double_t MCFM23NP_EYhstat_as[100];
     Double_t MCFM23NP_EYhstat_as_pdf[100];
     Double_t MCFM23NP_EYhstat_as_pdf_sca[100];
     Double_t MCFM23NP_EYhstat_as_pdf_sca_np[100];


     std::string hnameMCFM;
     std::string hnameaMCNLOmpi;
     std::string hnameaMCNLO;
     // graphical properties
     TLegend *legend = 0;  
     TLegend *legend1 = 0;  
     TLegend *legend2 = 0;  
     std::string xTitleStr;
     std::string yTitleStr;
     int colMCNLO = 2;
     int colfMCFMas  = 5;
     int colfMCFMas_pdf  = 2;
     int colfMCFMas_pdf_sca  = 94;
     int colfMCFMall = 6;
     int colSherpa   = 3;
     int colPythia   = 4;
     int collMCFM    = 1;
     int colmMCFM    = 1;
     int collMCFM8   = 216;
     int colmMCFM8   = 216;
     int collMCFM23  = 7;
     int colmMCFM23  = 7;
     
     std::string fdataname = "../data/Unfolded_"+hname+".txt";
     if (m_debug > 0) std::cout<<" fdataname will be "<<fdataname<<std::endl;
     

     const Int_t n0 = 2;
     Double_t xEmpty[n0], yEmpty[n0];
     for (Int_t i=0;i<n0;i++) {
       xEmpty[i] = (i+1)*0.00001;
       yEmpty[i] = 10000000.;
       //     printf(" i %i %f %f \n",i,x[i],y[i]);
     }
     TGraph* grEmpty = new TGraph(n0,xEmpty,yEmpty);
     

     TH2F* hempty; 

     Double_t dataNbbb[100];
     Double_t dataNBayes[100];
     Double_t dataMbbb[100];
     Double_t dataMBayes[100];
     Double_t edataNbbb[100];
     Double_t edataNBayes[100];
     Double_t edataMbbb[100];
     Double_t edataMBayes[100];
     for (int j=0; j<100; ++j) 
     {
	 dataNbbb[j] = 0.;
	 dataNBayes[j] = 0.;
	 dataMbbb[j] = 0.;
	 dataMBayes[j] = 0.;	 
	 edataNbbb[j] = 0.;
	 edataNBayes[j] = 0.;
	 edataMbbb[j] = 0.;
	 edataMBayes[j] = 0.;	 
     }
     


     if (hname=="j_bjet_Pt_binned")
       {
	 xTitleStr = "(b-jet selection) b-jet p_{T} [GeV]";
	 yTitleStr = "d#sigma/dp_{T} [fb/GeV]";
	 hnameMCFM = "Pt_gae_";
	 hnameaMCNLO = "bjet_pT";
	 if (oldaMCNLO) hnameaMCNLOmpi = "bjet_pT";
	 else hnameaMCNLOmpi = hname;
	 //grEmpty->GetYaxis()->SetRangeUser(0.01,500.);
	 if (setLogY) hempty = new TH2F("empty","empty",10,0.0,600.,10,.01,500.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hempty = new TH2F("empty","empty",10,0.0,600.,10,.0,300.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 //grEmpty->GetXaxis()->SetRangeUser(0.,600.);

	 

// 	 /*****.....**** Nicola unfolded yield  */
// 	 dataNbbb[0]=17889.2    ;
// 	 dataNbbb[1]=16641.7    ;
// 	 dataNbbb[2]=7705.4     ;
// 	 dataNbbb[3]=3704.6     ;
// 	 dataNbbb[4]=1720.5     ;
// 	 dataNbbb[5]=300.4      ;
// 	 ///////////////
// 	 /***********  Nicola rooUnfold */
//          dataNBayes[0]=17694.4  ;
//          dataNBayes[1]=16623	;
//          dataNBayes[2]=7726.9	;
//          dataNBayes[3]=3736.4	;
//          dataNBayes[4]=1705.7	;
//          dataNBayes[5]=301.8	;
// 	 ///////////////
// 	 /***** Mark bin by bin Alpgen */
// 	 dataMbbb[0]=	21230.8  ;
// 	 dataMbbb[1]=	16282.5	 ;
// 	 dataMbbb[2]=	7007.6	 ;
// 	 dataMbbb[3]=	3292	 ;
// 	 dataMbbb[4]=	1493.2	 ;
// 	 dataMbbb[5]=	259.3	 ;
// 	 ///////////////
// 	 /***** Mark bin by bin Sherpa*/
// 	 dataMBayes[0]=	21870.3  ;
// 	 dataMBayes[1]=	16553	 ;
// 	 dataMBayes[2]=	7299.7	 ;
// 	 dataMBayes[3]=	3488.9	 ;
// 	 dataMBayes[4]=	1675.4	 ;
// 	 dataMBayes[5]=	294.1	 ;
//	 ///////////////

       }
     else if ("j_bjet_Y_binned"==hname)
       {
	 xTitleStr = "(b-jet selection) b-jet |Y|";
	 yTitleStr = "d#sigma/d|Y| [fb/Y-unit]";
	 if (setLogY) hempty = new TH2F("empty","empty",10,0.0,3.,10,800.,4000.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hempty = new TH2F("empty","empty",10,0.0,3.,10,500.0,4000.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hnameMCFM = "Y_gae_";
	 hnameaMCNLO = "bjet_y";
	 if (oldaMCNLO) hnameaMCNLOmpi = "bjet_y";
	 else hnameaMCNLOmpi = hname;
       }
     else if ("j_bjet_ZY_binned"==hname)
       {
	 xTitleStr = "(Zb(b) selection) Z |Y|";
	 yTitleStr = "d#sigma/d|Y| [fb/Y-unit]";
	 if (setLogY) hempty = new TH2F("empty","empty",10,0.0,3.,10,100.,5000.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hempty = new TH2F("empty","empty",10,0.0,3.,10,.0,5000.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hnameMCFM = "ZY_gae_";
	 hnameaMCNLO = "Z_y";
	 if (oldaMCNLO) hnameaMCNLOmpi = "Z_y";
	 else hnameaMCNLOmpi = hname;
       }
     else if ("j_bjet_ZPt_binned"==hname)
       {
	 xTitleStr = "(Zb(b) selection) Z p_{T} [GeV]";
	 yTitleStr = "d#sigma/dp_{T} [fb/GeV]";
	 if (setLogY) hempty = new TH2F("empty","empty",10,0.0,600.,10,.03,300.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hempty = new TH2F("empty","empty",10,0.0,600.,10,.0,135.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hnameMCFM = "ZPt_gae_";
	 hnameaMCNLO = "Z_pT";
	 if (oldaMCNLO) hnameaMCNLOmpi = "Z_pT";
	 else hnameaMCNLOmpi = hname;
       }
     else if ("j_bjet_ZdR620_binned"==hname)
       {
	 xTitleStr = "(b-jet selection) #DeltaR(Z,b-jet)";
	 yTitleStr = "d#sigma/dR [fb/R-unit]";
	 if (setLogY) hempty = new TH2F("empty","empty",10,0.0,8.,10,.5,10000.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hempty = new TH2F("empty","empty",10,0.0,8.,10,.0,5500.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hnameMCFM = "ZPt_gae_";
	 hnameMCFM = "ZdR_gae_";
	 hnameaMCNLO = "Zj_dR";
	 if (oldaMCNLO) hnameaMCNLOmpi = "Zj_dR";
	 else hnameaMCNLOmpi = hname;	
       }
     else if ("j_bjet_ZdPhiF20_binned"==hname)
       {
	 xTitleStr = "(b-jet selection)  #Delta#phi(Z,b-jet)";
	 yTitleStr = "d#sigma/d#phi [fb/rad]";
	 if (setLogY) hempty = new TH2F("empty","empty",10,0.0,3.5,10,30.,20000.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hempty = new TH2F("empty","empty",10,0.0,3.5,10,.0,14000.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hnameMCFM = "ZdPhi_gae_";
	 hnameaMCNLO = "Zj_dphi";
	 if (oldaMCNLO) hnameaMCNLOmpi = "Zj_dphi";
	 else hnameaMCNLOmpi = hname;
       }
     else if ("j_2bjet_ZPt_binned"==hname)
       {
	 xTitleStr = "(Zbb selection) Z p_{T} [GeV]";
	 yTitleStr = "d#sigma/dp_{T} [fb/GeV]";
	 if (setLogY) hempty = new TH2F("empty","empty",10,0.0,300.,10,.2,10.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hempty = new TH2F("empty","empty",10,0.0,300.,10,.0,7.5); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hnameMCFM = "ZPT_gae_";
	 hnameaMCNLO = "Z_pT_bb";
	 if (oldaMCNLO) hnameaMCNLOmpi = "Z_pT_bb";
	 else hnameaMCNLOmpi = hname;
       }
     else if ("j_2bjet_dR_binned"==hname)
       {
	 xTitleStr = "(Zbb selection) #DeltaR(b-jet,b-jet)";
	 yTitleStr = "d#sigma/dR [fb/R-unit]";
	 if (setLogY) hempty = new TH2F("empty","empty",10,0.0,5.2,10,30.,300.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hempty = new TH2F("empty","empty",10,0.0,5.2,10,.0,300.); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hnameMCFM = "bbdR_gae_";
	 hnameaMCNLO = "bb_dR";
	 if (oldaMCNLO) hnameaMCNLOmpi = "bb_dR";
	 else hnameaMCNLOmpi = hname;
       }
     else if ("j_2bjet_Mbb_binned"==hname)
       {
	 xTitleStr = "(Zbb selection) M_{b-jet,b-jet}";
	 yTitleStr = "d#sigma/dM [fb/GeV]";
	 if (setLogY) hempty = new TH2F("empty","empty",10,0.0,400.,10,.2,7.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hempty = new TH2F("empty","empty",10,0.0,400.,10,.0,5.); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hnameMCFM = "bbM_gae_";
	 hnameaMCNLO = "bb_mass";
	 if (oldaMCNLO) hnameaMCNLOmpi = "bb_mass";
	 else hnameaMCNLOmpi = hname;
       }
     char* xTitle = xTitleStr.c_str();
     char* yTitle = yTitleStr.c_str();
     hempty->GetXaxis()->SetTitle(xTitle);
     hempty->GetYaxis()->SetTitle(yTitle);


     if (legend==0) {
       if ("j_2bjet_dR_binned"==hname || "j_bjet_ZdPhiF20_binned"==hname) legend=new TLegend(0.175,0.65,0.475,0.85);  
       else legend=new TLegend(0.6,0.65,0.95,0.85);  
     }
     legend->SetShadowColor(1);
     legend->SetBorderSize(0);
     legend->SetFillStyle(0);
     //legend->SetFillColor(kWhite);
     legend->SetTextSize(0.03);  
     if (legend1==0) {
       if ("j_2bjet_dR_binned"==hname || "j_bjet_ZdPhiF20_binned"==hname) legend1=new TLegend(0.175,0.55,0.475,0.90);  
       else legend1=new TLegend(0.53,0.55,0.98,0.90);  
     }
     legend1->SetShadowColor(1);
     legend1->SetBorderSize(0);
     legend1->SetFillStyle(0);
     //legend->SetFillColor(kWhite);
     legend1->SetTextSize(0.03);  
     if (legend2==0) {
       if ("j_2bjet_dR_binned"==hname || "j_bjet_ZdPhiF20_binned"==hname) legend2=new TLegend(0.175,0.55,0.475,0.90);  
       else legend2=new TLegend(0.53,0.55,0.98,0.90);  
     }
     legend2->SetShadowColor(1);
     legend2->SetBorderSize(0);
     legend2->SetFillStyle(0);
     //legend->SetFillColor(kWhite);
     legend2->SetTextSize(0.03);  
     
     // MCFM access 
     std::string hnameMCFMas;
     std::string hnameMCFMpdf;
     std::string hnameMCFMsta; 
     std::string hnameMCFMsca; 
     hnameMCFMas  = hnameMCFM+"as";
     hnameMCFMpdf = hnameMCFM+"pdf";
     hnameMCFMsta = hnameMCFM+"stat";
     hnameMCFMsca = hnameMCFM+"scale";
     // aMCatNLO access 

     if (m_debug > 0) std::cout<<"============ hnames MCFM/aMCNLO/aMCNLOmpi/Sherpa/Pythia= "<<hnameMCFMas<<"/"<<hnameaMCNLO<<"/"<<hnameaMCNLOmpi<<"/"<<hnameSherpa<<"/"<<hnamePythia<<std::endl;

     grMCFMas    = (TGraphAsymmErrors*)fileMCFM->Get(hnameMCFMas.c_str());
     grMCFMpdf   = (TGraphAsymmErrors*)fileMCFM->Get(hnameMCFMpdf.c_str());
     grMCFMsta   = (TGraphAsymmErrors*)fileMCFM->Get(hnameMCFMsta.c_str());
     grMCFMsca   = (TGraphAsymmErrors*)fileMCFM->Get(hnameMCFMsca.c_str());
     
     grMCFM8as    = (TGraphAsymmErrors*)fileMCFMmstw->Get(hnameMCFMas.c_str());
     grMCFM8pdf   = (TGraphAsymmErrors*)fileMCFMmstw->Get(hnameMCFMpdf.c_str());
     grMCFM8sta   = (TGraphAsymmErrors*)fileMCFMmstw->Get(hnameMCFMsta.c_str());
     grMCFM8sca   = (TGraphAsymmErrors*)fileMCFMmstw->Get(hnameMCFMsca.c_str());
     
     grMCFM23as    = (TGraphAsymmErrors*)fileMCFMnnpdf->Get(hnameMCFMas.c_str());
     grMCFM23pdf   = (TGraphAsymmErrors*)fileMCFMnnpdf->Get(hnameMCFMpdf.c_str());
     grMCFM23sta   = (TGraphAsymmErrors*)fileMCFMnnpdf->Get(hnameMCFMsta.c_str());
     grMCFM23sca   = (TGraphAsymmErrors*)fileMCFMnnpdf->Get(hnameMCFMsca.c_str());
     
     grMCFM   = grMCFMas;
     grMCFM8  = grMCFM8as;
     grMCFM23 = grMCFM23as;
     grSherpa = (TGraphAsymmErrors*)fileSHPY->Get(hnameSherpa.c_str());
     grPythia = (TGraphAsymmErrors*)fileSHPY->Get(hnamePythia.c_str());
     grSherpaNP = (TGraphAsymmErrors*)fileSHPY->Get(hnameSherpaNP.c_str());
     grPythiaNP = (TGraphAsymmErrors*)fileSHPY->Get(hnamePythiaNP.c_str());
     grNPSherpa = (TGraphAsymmErrors*)fileSHPY->Get(hnameNPSherpa.c_str());
     grNPPythia = (TGraphAsymmErrors*)fileSHPY->Get(hnameNPPythia.c_str());
     if (oldaMCNLO)
     {
         if (m_debug>2) std::cout<<" I believe I've asked for the old aMC@NLO ..... "<<std::endl;
	 TH1* haMCNLOMPI =  (TH1*)fileaMCNLOMPI->Get(hnameaMCNLOmpi.c_str());
	 if (hnameaMCNLO != "Z_y") graMCNLOMPI = new TGraphAsymmErrors(haMCNLOMPI);
     }
     else {
	 TH1F* haMCNLOMPI =  (TH1F*)fileaMCNLOMPI->Get(hnameaMCNLOmpi.c_str());
	 graMCNLOMPI =  new TGraphAsymmErrors(haMCNLOMPI); //(TGraphAsymmErrors*)fileaMCNLOMPI->Get(hnameaMCNLOmpi.c_str());
	 if (graMCNLOMPI!=NULL) if (m_debug>2) std::cout<<" got tgraph for amcnlo"<<std::endl;
	   else {
	     if (m_debug>2) std::cout<<" not found TGraph for amcnlo"<<std::endl;
	   }
     }     
     TH1* haMCNLO   =  (TH1*)fileaMCNLO->Get(hnameaMCNLO.c_str());
     if (hnameaMCNLO != "Z_y") graMCNLO   = new TGraphAsymmErrors(haMCNLO);
     if (m_debug>2) std::cout<< "histos/tgraph etc obtained also for aMC@NLO" <<std::endl;
     //return;
     

     //if (c1) delete c1;
     TCanvas* c1 = new TCanvas("NP","NP");
     grNPSherpa->SetLineColor(colSherpa);
     grNPSherpa->SetMarkerColor(colSherpa);
     grNPSherpa->Draw("APE");
     grNPPythia->SetLineColor(colPythia);
     grNPPythia->SetMarkerColor(colPythia);
     grNPPythia->Draw("PESAME");

     //return;
     
     Double_t* grNPSherpaY     = grNPSherpa->GetY();
     Double_t* grNPPythiaY     = grNPPythia->GetY();
     Double_t* grNPSherpaEYl   = grNPSherpa->GetEYlow();
     Double_t* grNPPythiaEYl   = grNPPythia->GetEYlow();
     Double_t* grNPSherpaEYh   = grNPSherpa->GetEYhigh();
     Double_t* grNPPythiaEYh   = grNPPythia->GetEYhigh();
     Double_t  grNPY[100]; 
     Double_t  grNPEYlstatAndSys[100]; 
     Double_t  grNPEYhstatAndSys[100]; 
     Double_t  grNPEYlstat[100]; 
     Double_t  grNPEYhstat[100]; 


     Double_t* grSherpaY       = grSherpa->GetY();
     Double_t* grSherpaEyLow   = grSherpa->GetEYlow();
     Double_t* grSherpaEyHigh  = grSherpa->GetEYhigh();
     Double_t  grSherpaYfb[100]; 
     Double_t  grSherpaEyLowfb[100];
     Double_t  grSherpaEyHighfb[100];
     Double_t  grSherpaNPYfb[100]; 
     Double_t  grSherpaNPEyLowfb[100];
     Double_t  grSherpaNPEyHighfb[100];
     Double_t* grPythiaY       = grPythia->GetY();
     Double_t* grPythiaEyLow   = grPythia->GetEYlow();
     Double_t* grPythiaEyHigh  = grPythia->GetEYhigh();
     Double_t  grPythiaYfb[100]; 
     Double_t  grPythiaEyLowfb[100];
     Double_t  grPythiaEyHighfb[100];
     Double_t  grPythiaNPYfb[100]; 
     Double_t  grPythiaNPEyLowfb[100];
     Double_t  grPythiaNPEyHighfb[100];
     Double_t* grSherpaX       = grSherpa->GetX();
     Double_t* grX             = grSherpa->GetX();
     Double_t* grEX            = grSherpa->GetEXlow();
     Double_t* graMCNLOMPIX;       if (hnameaMCNLO != "Z_y" || (!oldaMCNLO)) graMCNLOMPIX      = graMCNLOMPI->GetX();     
     Double_t* graMCNLOMPIY;       if (hnameaMCNLO != "Z_y" || (!oldaMCNLO)) graMCNLOMPIY      = graMCNLOMPI->GetY();
     Double_t* graMCNLOMPIEyLow;   if (hnameaMCNLO != "Z_y" || (!oldaMCNLO)) graMCNLOMPIEyLow  = graMCNLOMPI->GetEYlow();
     Double_t* graMCNLOMPIEyHigh;  if (hnameaMCNLO != "Z_y" || (!oldaMCNLO)) graMCNLOMPIEyHigh = graMCNLOMPI->GetEYhigh();
     Double_t* graMCNLOX;       if (hnameaMCNLO != "Z_y") graMCNLOX      = graMCNLO->GetX();
     Double_t* graMCNLOY;       if (hnameaMCNLO != "Z_y") graMCNLOY      = graMCNLO->GetY();
     Double_t* graMCNLOEyLow;   if (hnameaMCNLO != "Z_y") graMCNLOEyLow  = graMCNLO->GetEYlow();
     Double_t* graMCNLOEyHigh;  if (hnameaMCNLO != "Z_y") graMCNLOEyHigh = graMCNLO->GetEYhigh();
     Double_t  graMCNLOMPIYfb[100]; 
     Double_t  graMCNLOMPIEyLowfb[100];
     Double_t  graMCNLOMPIEyHighfb[100];
     Double_t  graMCNLOYfb[100]; 
     Double_t  graMCNLOEyLowfb[100];
     Double_t  graMCNLOEyHighfb[100];

     //return;
     if (m_debug>1) std::cout<<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% reading data"<<std::endl;     
     Double_t X1[100], X2[100], X3[100], X4[100];
     Int_t nrows = 0;
     if (hname.substr(0,6)=="j_bjet") {
       readVectorsFromFile(fdataname,nrows,X1,X2,X3,X4,m_debug);
       for (int j=0; j<nrows; ++j)
	 {
	   dataNbbb[j]   = X1[j];
	   dataNBayes[j] = X2[j];
	   dataMbbb[j]   = X3[j];
	   dataMBayes[j] = X4[j];
	   if (m_debug>1) if (m_debug>1) std::cout<<" data =="<<dataNBayes[j] <<std::endl;
	 }
     }
     //return;
     if (m_debug>1) std::cout<<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% data read"<<std::endl;

     
     int Npoints = grSherpa->GetN();
     int offset = 0;
     if (hname=="j_bjet_Pt_binned" && oldaMCNLO) offset=2;
     for (int k=0; k<Npoints; ++k)
     {
	 Double_t dx = 2.*grEX[k];

	 if (hname.substr(0,6)=="j_bjet") {
	   // data 
	   double aa = dataNbbb[k];
	   edataNbbb[k]=sqrt(dataNbbb[k]);
	   edataMbbb[k]=sqrt(dataMbbb[k]);
	   edataNBayes[k]=sqrt(dataNBayes[k]);
	   edataMBayes[k]=sqrt(dataMBayes[k]);
	   
	   dataNbbb[k]=  0.5 * dataNbbb[k]  /(dx*Lumi);
	   dataMbbb[k]=  0.5 * dataMbbb[k]  /(dx*Lumi);
	   dataNBayes[k]=0.5 * dataNBayes[k]/(dx*Lumi);
	   dataMBayes[k]=0.5 * dataMBayes[k]/(dx*Lumi);
	   //
	   edataNbbb[k]  =0.5 * edataNbbb[k]/(dx*Lumi);
	   edataMbbb[k]  =0.5 * edataMbbb[k]/(dx*Lumi);
	   edataNBayes[k]=0.5 * edataNBayes[k]/(dx*Lumi);
	   edataMBayes[k]=0.5 * edataMBayes[k]/(dx*Lumi);
	   if (m_debug>1) std::cout<<"Data ....... "<<aa<<" after scaling in bin 1 "<<dataNbbb[k]<<" dx = "<<dx<<" sherpa "<<grNPSherpaY[k]<<std::endl;
	 }
	 // NP 
	 grNPY[k]  = 0.5*(grNPSherpaY[k]+grNPPythiaY[k]); // multiplicative term
	 grNPEYlstat[k]  = sqrt(grNPSherpaEYl[k]*grNPSherpaEYl[k]+grNPPythiaEYl[k]*grNPPythiaEYl[k]); // multiplicative term
	 grNPEYhstat[k]  = sqrt(grNPSherpaEYh[k]*grNPSherpaEYh[k]+grNPPythiaEYh[k]*grNPPythiaEYh[k]); // multiplicative term
	 double grNPEYsys = 0.5*abs(grNPSherpaY[k]-grNPPythiaY[k]); // multiplicative term
	 grNPEYlstatAndSys[k]  = sqrt(grNPEYlstat[k]*grNPEYlstat[k] + grNPEYsys*grNPEYsys); // multiplicative term
	 grNPEYhstatAndSys[k]  = sqrt(grNPEYhstat[k]*grNPEYhstat[k] + grNPEYsys*grNPEYsys); // multiplicative term
	 if (m_debug>0) {
	   std::cout<<"NPtest::::::  SH, PY, statSH, statPY "<<grNPSherpaY[k]<<" "<<grNPPythiaY[k]<<" "<<grNPSherpaEYl[k]<<" "<<grNPPythiaEYl[k]<<std::endl;
	   std::cout<<"NPtest::::::  NP, NP_estat, NP_esyst, NP_etot "<<grNPY[k]<<" "<<grNPEYlstat[k]<<" "<<grNPEYsys<<" "<<grNPEYhstatAndSys[k]<<std::endl;
	 }

	 //MCFM corrected 
         MCFMNP_Y[k]       = grNPY[k] * grMCFM->GetY()[k];
         MCFM8NP_Y[k]      = grNPY[k] * grMCFM8->GetY()[k];
         MCFM23NP_Y[k]     = grNPY[k] * grMCFM23->GetY()[k];
	 
	 if (m_debug>0) {
	   std::cout<<"NP and MCFM ct10     test :::::::::: beforeNP, after NP "<<grMCFM->GetY()[k]  <<" "<<MCFMNP_Y[k]<<std::endl;
	   std::cout<<"NP and MCFM mstw2008 test :::::::::: beforeNP, after NP "<<grMCFM8->GetY()[k] <<" "<<MCFM8NP_Y[k]<<std::endl;
	   std::cout<<"NP and MCFM nnpdf23  test :::::::::: beforeNP, after NP "<<grMCFM23->GetY()[k]<<" "<<MCFM23NP_Y[k]<<std::endl;
	 }
	 MCFMNP_Y[k]       = 0.97* MCFMNP_Y[k];  // after FSR 
	 MCFM8NP_Y[k]      = 0.97* MCFM8NP_Y[k]; // after FSR 
	 MCFM23NP_Y[k]     = 0.97* MCFM23NP_Y[k];// after FSR 
	 
	 double staErrL    = grMCFMsta->GetEYlow()[k];
	 double staErrH    = grMCFMsta->GetEYhigh()[k];
	 double asErrL  = grMCFMas ->GetEYlow()[k];
	 double pdfErrL = grMCFMpdf->GetEYlow()[k];
	 double scaErrL = grMCFMsca->GetEYlow()[k];
	 double asErrH  = grMCFMas ->GetEYhigh()[k];
	 double pdfErrH = grMCFMpdf->GetEYhigh()[k];
	 double scaErrH = grMCFMsca->GetEYhigh()[k];
	 double npErrL  = grNPEYlstatAndSys[k]*MCFMNP_Y[k];// multiplicative term
	 double npErrH  = grNPEYhstatAndSys[k]*MCFMNP_Y[k];// multiplicative term
	 npErrL  = sqrt(npErrL*npErrL + MCFMNP_Y[k]*MCFMNP_Y[k]*0.0001); // multiplicative errors from NP and 1% for FSR
	 npErrH  = sqrt(npErrH*npErrH + MCFMNP_Y[k]*MCFMNP_Y[k]*0.0001); // multiplicative errors from NP and 1% for FSR
	 if (m_debug>0) std::cout<<"NP and MCFM test :::::::::: npErrL, H  "<<npErrL<<" "<<npErrH<<std::endl;

         MCFMNP_EYlstat[k] = staErrL;
         MCFMNP_EYlstat_as[k]            = sqrt(MCFMNP_EYlstat[k]       *MCFMNP_EYlstat[k]                + 
						asErrL*asErrL);
         MCFMNP_EYlstat_as_pdf[k]        = sqrt(MCFMNP_EYlstat_as[k]    *MCFMNP_EYlstat_as[k]             + 
						pdfErrL*pdfErrL);
         MCFMNP_EYlstat_as_pdf_sca[k]    = sqrt(MCFMNP_EYlstat_as_pdf[k]*MCFMNP_EYlstat_as_pdf[k]         + 
						scaErrL*scaErrL);
         MCFMNP_EYlstat_as_pdf_sca_np[k] = sqrt(MCFMNP_EYlstat_as_pdf_sca[k]*MCFMNP_EYlstat_as_pdf_sca[k] + 
						npErrL*npErrL);
         MCFMNP_EYhstat[k] = staErrH;
         MCFMNP_EYhstat_as[k]            = sqrt(MCFMNP_EYhstat[k]       *MCFMNP_EYhstat[k]                + 
						asErrH *asErrH);
         MCFMNP_EYhstat_as_pdf[k]        = sqrt(MCFMNP_EYhstat_as[k]    *MCFMNP_EYhstat_as[k]             + 
						pdfErrH*pdfErrH);
         MCFMNP_EYhstat_as_pdf_sca[k]    = sqrt(MCFMNP_EYhstat_as_pdf[k]*MCFMNP_EYhstat_as_pdf[k]         + 
						scaErrH*scaErrH);
         MCFMNP_EYhstat_as_pdf_sca_np[k] = sqrt(MCFMNP_EYhstat_as_pdf_sca[k]*MCFMNP_EYhstat_as_pdf_sca[k] + 
						npErrH*npErrH); 
	 if (m_debug>0) {
	   std::cout<<"MCFMwNP errors test :::::::::: stat, a_s, pdf, sca, np (low)  "<<staErrL<<" "<<asErrL<<" "<<pdfErrL<<" "<<scaErrL<<" "<<npErrL<<" tot "<<MCFMNP_EYlstat_as_pdf_sca_np[k]<<std::endl;
	   std::cout<<"MCFMwNP errors test :::::::::: stat, a_s, pdf, sca, np (high) "<<staErrH<<" "<<asErrH<<" "<<pdfErrH<<" "<<scaErrH<<" "<<npErrH<<" tot "<<MCFMNP_EYhstat_as_pdf_sca_np[k]<<std::endl;
	 }

	 staErrL    = sqrt(5.)*staErrL;
	 staErrH    = sqrt(5.)*staErrH;
	 asErrL  = grMCFM8as ->GetEYlow()[k];
	 pdfErrL = grMCFM8pdf->GetEYlow()[k];
	 scaErrL = grMCFM8sca->GetEYlow()[k];
	 asErrH  = grMCFM8as ->GetEYhigh()[k];
	 pdfErrH = grMCFM8pdf->GetEYhigh()[k];
	 scaErrH = grMCFM8sca->GetEYhigh()[k];
	 npErrL  = grNPEYlstatAndSys[k]*MCFM8NP_Y[k];// multiplicative term
	 npErrH  = grNPEYhstatAndSys[k]*MCFM8NP_Y[k];// multiplicative term
	 npErrL  = sqrt(npErrL*npErrL + MCFM8NP_Y[k]*MCFM8NP_Y[k]*0.0001); // multiplicative errors from NP and 1% for FSR
	 npErrH  = sqrt(npErrH*npErrH + MCFM8NP_Y[k]*MCFM8NP_Y[k]*0.0001); // multiplicative errors from NP and 1% for FSR

	 if (m_debug>0) std::cout<<"NP and MCFM test :::::::::: npErrL, H  "<<npErrL<<" "<<npErrH<<std::endl;

         MCFM8NP_EYlstat[k] = staErrL;
         MCFM8_EYlstat[k]   = staErrL;
         MCFM8NP_EYlstat_as[k]            = sqrt(MCFM8NP_EYlstat[k]       *MCFM8NP_EYlstat[k]                + 
						 asErrL*asErrL);	
         MCFM8NP_EYlstat_as_pdf[k]        = sqrt(MCFM8NP_EYlstat_as[k]    *MCFM8NP_EYlstat_as[k]             + 
						pdfErrL*pdfErrL);
         MCFM8NP_EYlstat_as_pdf_sca[k]    = sqrt(MCFM8NP_EYlstat_as_pdf[k]*MCFM8NP_EYlstat_as_pdf[k]         + 
						scaErrL*scaErrL);
         MCFM8NP_EYlstat_as_pdf_sca_np[k] = sqrt(MCFM8NP_EYlstat_as_pdf_sca[k]*MCFM8NP_EYlstat_as_pdf_sca[k] + 
						npErrL*npErrL);
         MCFM8NP_EYhstat[k] = staErrH;
         MCFM8_EYhstat[k]   = staErrH;
         MCFM8NP_EYhstat_as[k]            = sqrt(MCFM8NP_EYhstat[k]       *MCFM8NP_EYhstat[k]                + 
						asErrH *asErrH);
         MCFM8NP_EYhstat_as_pdf[k]        = sqrt(MCFM8NP_EYhstat_as[k]    *MCFM8NP_EYhstat_as[k]             + 
						pdfErrH*pdfErrH);
         MCFM8NP_EYhstat_as_pdf_sca[k]    = sqrt(MCFM8NP_EYhstat_as_pdf[k]*MCFM8NP_EYhstat_as_pdf[k]         + 
						scaErrH*scaErrH);
         MCFM8NP_EYhstat_as_pdf_sca_np[k] = sqrt(MCFM8NP_EYhstat_as_pdf_sca[k]*MCFM8NP_EYhstat_as_pdf_sca[k] + 
						npErrH*npErrH); 
	 if (m_debug>0) {
	   std::cout<<"MCFM8wNP errors test :::::::::: stat, a_s, pdf, sca, np (low)  "<<staErrL<<" "<<asErrL<<" "<<pdfErrL<<" "<<scaErrL<<" "<<npErrL<<" tot "<<MCFM8NP_EYlstat_as_pdf_sca_np[k]<<std::endl;
	   std::cout<<"MCFM8wNP errors test :::::::::: stat, a_s, pdf, sca, np (high) "<<staErrH<<" "<<asErrH<<" "<<pdfErrH<<" "<<scaErrH<<" "<<npErrH<<" tot "<<MCFM8NP_EYhstat_as_pdf_sca_np[k]<<std::endl;
	 }

	 //	 staErrL    = sqrt(5.)*staErrL;
	 //	 staErrH    = sqrt(5.)*staErrH;
	 asErrL  = grMCFM23as ->GetEYlow()[k];
	 pdfErrL = grMCFM23pdf->GetEYlow()[k];
	 scaErrL = grMCFM23sca->GetEYlow()[k];
	 asErrH  = grMCFM23as ->GetEYhigh()[k];
	 pdfErrH = grMCFM23pdf->GetEYhigh()[k];
	 scaErrH = grMCFM23sca->GetEYhigh()[k];
	 npErrL  = grNPEYlstatAndSys[k]*MCFM23NP_Y[k];// multiplicative term
	 npErrH  = grNPEYhstatAndSys[k]*MCFM23NP_Y[k];// multiplicative term
	 npErrL  = sqrt(npErrL*npErrL + MCFM23NP_Y[k]*MCFM23NP_Y[k]*0.0001); // multiplicative errors from NP and 1% for FSR
	 npErrH  = sqrt(npErrH*npErrH + MCFM23NP_Y[k]*MCFM23NP_Y[k]*0.0001); // multiplicative errors from NP and 1% for FSR

	 if (m_debug>0) std::cout<<"NP and MCFM test :::::::::: npErrL, H  "<<npErrL<<" "<<npErrH<<std::endl;

         MCFM23NP_EYlstat[k] = MCFM8NP_EYlstat[k];
         MCFM23_EYlstat[k]   = MCFM8_EYlstat[k];
         MCFM23NP_EYlstat_as[k]            = sqrt(MCFM23NP_EYlstat[k]       *MCFM23NP_EYlstat[k]                + 
						 asErrL*asErrL);	
         MCFM23NP_EYlstat_as_pdf[k]        = sqrt(MCFM23NP_EYlstat_as[k]    *MCFM23NP_EYlstat_as[k]             + 
						pdfErrL*pdfErrL);
         MCFM23NP_EYlstat_as_pdf_sca[k]    = sqrt(MCFM23NP_EYlstat_as_pdf[k]*MCFM23NP_EYlstat_as_pdf[k]         + 
						scaErrL*scaErrL);
         MCFM23NP_EYlstat_as_pdf_sca_np[k] = sqrt(MCFM23NP_EYlstat_as_pdf_sca[k]*MCFM23NP_EYlstat_as_pdf_sca[k] + 
						npErrL*npErrL);
         MCFM23NP_EYhstat[k] = MCFM8NP_EYhstat[k];
         MCFM23_EYhstat[k]   = MCFM8_EYhstat[k];
         MCFM23NP_EYhstat_as[k]            = sqrt(MCFM23NP_EYhstat[k]       *MCFM23NP_EYhstat[k]                + 
						asErrH *asErrH);
         MCFM23NP_EYhstat_as_pdf[k]        = sqrt(MCFM23NP_EYhstat_as[k]    *MCFM23NP_EYhstat_as[k]             + 
						pdfErrH*pdfErrH);
         MCFM23NP_EYhstat_as_pdf_sca[k]    = sqrt(MCFM23NP_EYhstat_as_pdf[k]*MCFM23NP_EYhstat_as_pdf[k]         + 
						scaErrH*scaErrH);
         MCFM23NP_EYhstat_as_pdf_sca_np[k] = sqrt(MCFM23NP_EYhstat_as_pdf_sca[k]*MCFM23NP_EYhstat_as_pdf_sca[k] + 
						npErrH*npErrH); 
	 if (m_debug>0) {
	   std::cout<<"MCFM23wNP errors test :::::::::: stat, a_s, pdf, sca, np (low)  "<<staErrL<<" "<<asErrL<<" "<<pdfErrL<<" "<<scaErrL<<" "<<npErrL<<" tot "<<MCFM23NP_EYlstat_as_pdf_sca_np[k]<<std::endl;
	   std::cout<<"MCFM23wNP errors test :::::::::: stat, a_s, pdf, sca, np (high) "<<staErrH<<" "<<asErrH<<" "<<pdfErrH<<" "<<scaErrH<<" "<<npErrH<<" tot "<<MCFM23NP_EYhstat_as_pdf_sca_np[k]<<std::endl;
	 }


	 grSherpaYfb[k]      = 1000.*grSherpaY[k];
	 grSherpaEyLowfb[k]  = 1000.*grSherpaEyLow[k];
	 grSherpaEyHighfb[k] = 1000.*grSherpaEyHigh[k];
	 grPythiaYfb[k]      = 1000.*grPythiaY[k];
	 grPythiaEyLowfb[k]  = 1000.*grPythiaEyLow[k];
	 grPythiaEyHighfb[k] = 1000.*grPythiaEyHigh[k];

	 grSherpaNPYfb[k]      = 1000.*grSherpaNP->GetY()[k];	  
	 grSherpaNPEyLowfb[k]  = 1000.*grSherpaNP->GetEYlow()[k]; 
	 grSherpaNPEyHighfb[k] = 1000.*grSherpaNP->GetEYhigh()[k];
	 grPythiaNPYfb[k]      = 1000.*grPythiaNP->GetY()[k];	  
	 grPythiaNPEyLowfb[k]  = 1000.*grPythiaNP->GetEYlow()[k]; 
	 grPythiaNPEyHighfb[k] = 1000.*grPythiaNP->GetEYhigh()[k];

	 if (hnameaMCNLO != "Z_y" && oldaMCNLO)
	 {
	     graMCNLOYfb[k]      = 1000000.*graMCNLOY[k+offset];
	     graMCNLOEyLowfb[k]  = 1000000.*graMCNLOEyLow[k+offset];
	     graMCNLOEyHighfb[k] = 1000000.*graMCNLOEyHigh[k+offset];
	     graMCNLOMPIYfb[k]      = 1000.*graMCNLOMPIY[k+offset];
	     graMCNLOMPIEyLowfb[k]  = 1000.*graMCNLOMPIEyLow[k+offset];
	     graMCNLOMPIEyHighfb[k] = 1000.*graMCNLOMPIEyHigh[k+offset];
	 }
	 else if (!oldaMCNLO)
	 {
	     graMCNLOMPIYfb[k]      = mpiGlobCorr*100.*graMCNLOMPIY[k+offset];     ////// why 100. ???? - ok 1.0286232 is ZPt wMPI/woMPI total
	     graMCNLOMPIEyLowfb[k]  = mpiGlobCorr*100.*graMCNLOMPIEyLow[k+offset]; ////// why 100. ???? - ok 
	     graMCNLOMPIEyHighfb[k] = mpiGlobCorr*100.*graMCNLOMPIEyHigh[k+offset];////// why 100. ???? - ok 
	     if (hnameaMCNLO != "Z_y")
	     {
		 graMCNLOYfb[k]      = 1000000.*graMCNLOY[k+offset];
		 graMCNLOEyLowfb[k]  = 1000000.*graMCNLOEyLow[k+offset];
		 graMCNLOEyHighfb[k] = 1000000.*graMCNLOEyHigh[k+offset];
	     }	
	 }
     }
     TGraphAsymmErrors *grdataNbbb;
     if (hname.substr(0,6)=="j_bjet") grdataNbbb  = new TGraphAsymmErrors(Npoints, 
									  grX,                             dataNbbb, 
									  grEX,                            grEX,
									  edataNbbb,                       edataNbbb);
     TGraphAsymmErrors *grdataMbbb;
     if (hname.substr(0,6)=="j_bjet") grdataMbbb  = new TGraphAsymmErrors(Npoints, 
									  grX,                             dataMbbb, 
									  grEX,                            grEX,
									  edataMbbb,                       edataMbbb);
     TGraphAsymmErrors *grdataNBayes;
     if (hname.substr(0,6)=="j_bjet") grdataNBayes = new TGraphAsymmErrors(Npoints, 
									   grX,                             dataNBayes, 
									   grEX,                            grEX,
									   edataNBayes,                     edataNBayes);
     TGraphAsymmErrors *grdataMBayes;
     if (hname.substr(0,6)=="j_bjet") grdataMBayes = new TGraphAsymmErrors(Npoints, 
									   grX,                             dataMBayes, 
									   grEX,                            grEX,
									   edataMBayes,                     edataMBayes);
     TGraphAsymmErrors *grMCFMsta  = new TGraphAsymmErrors(Npoints, 
							   grX,                             grMCFM->GetY(), 
							   grEX,                            grEX,
							   grMCFM->GetEYlow(),              grMCFM->GetEYhigh());
     TGraphAsymmErrors *grMCFM8sta = new TGraphAsymmErrors(Npoints, 
							   grX,                             grMCFM8->GetY(), 
							   grEX,                            grEX,
							   MCFM8_EYlstat,                   MCFM8_EYhstat);
     TGraphAsymmErrors *grMCFM23sta = new TGraphAsymmErrors(Npoints, 
							    grX,                             grMCFM23->GetY(), 
							    grEX,                            grEX,
							    MCFM23_EYlstat,                  MCFM23_EYhstat);
     grMCFMNPsta               = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFMNP_Y, 
						       grEX,                            grEX,
						       MCFMNP_EYlstat,                  MCFMNP_EYhstat);
     grMCFM8NPsta              = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFM8NP_Y, 
						       grEX,                            grEX,
						       MCFM8NP_EYlstat,                 MCFM8NP_EYhstat);
     grMCFM23NPsta             = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFM23NP_Y, 
						       grEX,                            grEX,
						       MCFM23NP_EYlstat,                MCFM23NP_EYhstat);
     grMCFMNPsta_as            = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFMNP_Y, 
						       grEX,                            grEX,
						       MCFMNP_EYlstat_as,               MCFMNP_EYhstat_as);
     grMCFM8NPsta_as           = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFM8NP_Y, 
						       grEX,                            grEX,
						       MCFM8NP_EYlstat_as,              MCFM8NP_EYhstat_as);
     grMCFM23NPsta_as          = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFM23NP_Y, 
						       grEX,                            grEX,
						       MCFM23NP_EYlstat_as,             MCFM23NP_EYhstat_as);
     grMCFMNPsta_as_pdf        = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFMNP_Y, 
						       grEX,                            grEX,
						       MCFMNP_EYlstat_as_pdf,           MCFMNP_EYhstat_as_pdf);
     grMCFM8NPsta_as_pdf       = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFM8NP_Y, 
						       grEX,                            grEX,
						       MCFM8NP_EYlstat_as_pdf,          MCFM8NP_EYhstat_as_pdf);
     grMCFM23NPsta_as_pdf      = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFM23NP_Y, 
						       grEX,                            grEX,
						       MCFM23NP_EYlstat_as_pdf,         MCFM23NP_EYhstat_as_pdf);
     grMCFMNPsta_as_pdf_sca    = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFMNP_Y, 
						       grEX,                            grEX,
						       MCFMNP_EYlstat_as_pdf_sca,       MCFMNP_EYhstat_as_pdf_sca);
     grMCFM8NPsta_as_pdf_sca   = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFM8NP_Y, 
						       grEX,                            grEX,
						       MCFM8NP_EYlstat_as_pdf_sca,      MCFM8NP_EYhstat_as_pdf_sca);
     grMCFM23NPsta_as_pdf_sca  = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFM23NP_Y, 
						       grEX,                            grEX,
						       MCFM23NP_EYlstat_as_pdf_sca,     MCFM23NP_EYhstat_as_pdf_sca);
     grMCFMNPsta_as_pdf_sca_np = new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFMNP_Y, 
						       grEX,                            grEX,
						       MCFMNP_EYlstat_as_pdf_sca_np,    MCFMNP_EYhstat_as_pdf_sca_np);
     grMCFM8NPsta_as_pdf_sca_np= new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFM8NP_Y, 
						       grEX,                            grEX,
						       MCFM8NP_EYlstat_as_pdf_sca_np,   MCFM8NP_EYhstat_as_pdf_sca_np);
     grMCFM23NPsta_as_pdf_sca_np= new TGraphAsymmErrors(Npoints, 
						       grX,                             MCFM23NP_Y, 
						       grEX,                            grEX,
						       MCFM23NP_EYlstat_as_pdf_sca_np,  MCFM23NP_EYhstat_as_pdf_sca_np);

     TGraphAsymmErrors *grSherpafb = new TGraphAsymmErrors(Npoints, 
							   grX,                  grSherpaYfb, 
							   grEX,                 grEX,
							   grSherpaEyLowfb,      grSherpaEyHighfb);
     TGraphAsymmErrors *grSherpaNPfb = new TGraphAsymmErrors(Npoints, 
							     grX,                  grSherpaNPYfb, 
							     grEX,                 grEX,
							     grSherpaNPEyLowfb,    grSherpaNPEyHighfb);
     TGraphAsymmErrors *grPythiafb = new TGraphAsymmErrors(Npoints, 
							   grPythia->GetX(),     grPythiaYfb, 
							   grEX,                 grEX,
							   grPythiaEyLowfb,      grPythiaEyHighfb);
     TGraphAsymmErrors *grPythiaNPfb = new TGraphAsymmErrors(Npoints, 
							     grPythia->GetX(),     grPythiaNPYfb, 
							     grEX,                 grEX,
							     grPythiaNPEyLowfb,    grPythiaNPEyHighfb);

     std::cout<<"=============================================="<<std::endl;
     std::cout<<"====== Xsec Summary = "<<hname<<std::endl;
     std::cout<<"=============================================="<<std::endl;
     std::cout<<"====== MCFM     CT10 sta err NO NP ==========="<<std::endl;
     Double_t xsMCFMsta                                            = ComputeTotalCrossSection(grMCFMsta);
     std::cout<<"====== MCFM MSTW2008 sta err NO NP ==========="<<std::endl;
     Double_t xsMCFM8sta                                            = ComputeTotalCrossSection(grMCFM8sta);
     std::cout<<"====== MCFM NNPDF23  sta err NO NP ==========="<<std::endl;
     Double_t xsMCFM23sta                                           = ComputeTotalCrossSection(grMCFM23sta);
     std::cout<<"=============================================="<<std::endl;
     std::cout<<"====== SHERPA no NP =========================="<<std::endl;
     Double_t xsSHERPA                                             = ComputeTotalCrossSection(grSherpafb);
     std::cout<<"=============================================="<<std::endl;
     std::cout<<"====== PYTHIA no NP =========================="<<std::endl;
     Double_t xsPYTHIA                                             = ComputeTotalCrossSection(grPythiafb);
     std::cout<<"=============================================="<<std::endl;
     std::cout<<"====== MCFM sta err =========================="<<std::endl;
     Double_t xsMCFMNPsta                                          = ComputeTotalCrossSection(grMCFMNPsta);
     std::cout<<"====== MCFM sta_as err ======================="<<std::endl;
     Double_t xsMCFMNPsta_as                                       = ComputeTotalCrossSection(grMCFMNPsta_as);
     std::cout<<"====== MCFM sta_as_pdf err ==================="<<std::endl;
     Double_t xsMCFMNPsta_as_pdf                                   = ComputeTotalCrossSection(grMCFMNPsta_as_pdf);
     std::cout<<"====== MCFM sta_as_pdf_sca err ==============="<<std::endl;
     Double_t xsMCFMNPsta_as_pdf_sca                               = ComputeTotalCrossSection(grMCFMNPsta_as_pdf_sca);
     std::cout<<"====== MCFM CT10     sta_as_pdf_sca_np err ==="<<std::endl;
     Double_t xsMCFMNPsta_as_pdf_sca_np                          = ComputeTotalCrossSection(grMCFMNPsta_as_pdf_sca_np);
     std::cout<<"====== MCFM MSTW2008 sta_as_pdf_sca_np err ==="<<std::endl;
     Double_t xsMCFM8NPsta_as_pdf_sca_np                         = ComputeTotalCrossSection(grMCFM8NPsta_as_pdf_sca_np);
     std::cout<<"====== MCFM NNPDF23  sta_as_pdf_sca_np err ==="<<std::endl;
     Double_t xsMCFM23NPsta_as_pdf_sca_np                        = ComputeTotalCrossSection(grMCFM23NPsta_as_pdf_sca_np);

     Double_t xsaMCNLOMPI = 0.;
     Double_t xsaMCNLO   = 0.;
     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO)) {
       graMCNLOMPIfb = new TGraphAsymmErrors(grSherpa->GetN(), 
					    grX,                  graMCNLOMPIYfb, 
					    grEX,                 grEX,
					    graMCNLOMPIEyLowfb,    graMCNLOMPIEyHighfb);
       graMCNLOfb   = new TGraphAsymmErrors(grSherpa->GetN(), 
					    grX,                  graMCNLOYfb, 
					    grEX,                 grEX,
					    graMCNLOEyLowfb,      graMCNLOEyHighfb);
       std::cout<<"=============================================="<<std::endl; 
       std::cout<<"====== aMC@NLO no MPI ========================"<<std::endl;
       xsaMCNLO                                                     = ComputeTotalCrossSection(graMCNLOfb);
       std::cout<<"====== aMC@NLO NP ============================"<<std::endl;
       xsaMCNLOMPI                                                   = ComputeTotalCrossSection(graMCNLOMPIfb);
     }
     std::cout<<"=============================================="<<std::endl;
     std::cout<<"====== SHERPA NP ============================="<<std::endl;
     Double_t xsSHERPANP                                           = ComputeTotalCrossSection(grSherpaNPfb);
     std::cout<<"=============================================="<<std::endl;
     std::cout<<"====== PYTHIA NP ============================="<<std::endl;
     Double_t xsPYTHIANP                                           = ComputeTotalCrossSection(grPythiaNPfb);
     std::cout<<"=============================================="<<std::endl;
     std::cout<<"=============================================="<<std::endl;
     if (hname.substr(0,6)=="j_bjet" && writeDataResults) 
       {
	 std::cout<<"=============================================="<<std::endl;
	 std::cout<<"====== NIKO bxb =============================="<<std::endl;
	 Double_t xsDataNbxb   = ComputeTotalCrossSection(grdataNbbb);
	 std::cout<<"=============================================="<<std::endl;
	 std::cout<<"====== NIKO Bayes ============================"<<std::endl;
	 Double_t xsDataNBaye  = ComputeTotalCrossSection(grdataNBayes);
	 std::cout<<"=============================================="<<std::endl;
	 std::cout<<"====== MARK bxb =============================="<<std::endl;
	 Double_t xsDataMbxb   = ComputeTotalCrossSection(grdataMbbb);
	 std::cout<<"=============================================="<<std::endl;
	 std::cout<<"====== MARK Bayes ============================"<<std::endl;
	 Double_t xsDataMBayes = ComputeTotalCrossSection(grdataMBayes);
	 std::cout<<"=============================================="<<std::endl;
	 std::cout<<"=============================================="<<std::endl;
       } 

     //return;
     

     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO)) {
	 //if (cPlotNP2) delete cPlotNP2;
	 
     TCanvas *cPlotNP2 = new TCanvas("aMCNLOcheck","aMCNLO check");

     cPlotNP2->Divide(1,2);
     cPlotNP2->cd(1);
     graMCNLOMPI->SetLineColor(colMCNLO);
     graMCNLOMPI->SetMarkerColor(colMCNLO);
     graMCNLOMPI->Draw("APC");
     cPlotNP2->cd(2);
     if (graMCNLO) 
     {
	 graMCNLO->SetLineColor(colMCNLO);
	 graMCNLO->SetMarkerColor(colMCNLO);
	 graMCNLO->Draw("APC");
     }
     if (m_debug>1){
     for (int j=0; j<graMCNLOMPI->GetN(); ++j)
       {
	 std::cout<<"------- aMC@NLO orig TGRAPH j="<<j<<" x, y "<<graMCNLOMPIfb->GetX()[j]<<" "<<graMCNLOMPIfb->GetY()[j]<<" "<<graMCNLOfb->GetY()[j]<<std::endl;
       }
     }
     }
     std::string fname;
     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))  
       {
	   // if (cPlot2) delete cPlot2;
	   TCanvas *cPlot2 = new TCanvas("aMC@LNO","aMC@LNO");
	 
	 graMCNLOMPIfb->SetLineColor(2);
	 graMCNLOMPIfb->SetMarkerColor(2);
	 graMCNLOMPIfb->Draw("APC");
	 if (showMCNLOnoMPI) 
	 {
	     graMCNLOfb->SetLineColor(2);
	     graMCNLOfb->SetLineStyle(2);
	     graMCNLOfb->SetMarkerColor(2);
	     graMCNLOfb->SetMarkerStyle(24);
	     graMCNLOfb->Draw("PCSAME");
	 }
	 
	 legend->Clear();
	 legend->AddEntry(graMCNLOMPIfb, "aMC@NLO","PEL");	
	 if (showMCNLOnoMPI) legend->AddEntry(graMCNLOfb,   "aMC@NLO no MPI","PEL");	
	 legend->Draw();
	 if (saveAsEps)
	   {
	     fname = "fig"+hname+"_MCNLO.eps";
	     if (setLogY) fname = "log"+fname;
	     cPlot2->SaveAs(fname.c_str());
	   }
	 // if (saveAsRoot)
	 //   {
	 //     fname = "fig"+hname+"_MCNLO.root";
	 //     cPlot2->SaveAs(fname.c_str());
	 //   }
       }

     //return;
     //if (cPlot) delete cPlot;
     TCanvas *cPlot = new TCanvas("MCFM_SherpaPythia","MCFM_SherpaPythia");
     if (setLogY) 
       {
	 cPlot->SetLogy();
       }


     hempty->Draw();

     grMCFMsta->SetLineColor(collMCFM);
     grMCFMsta->SetMarkerColor(colmMCFM);
     grMCFMsta->SetLineWidth(1);
     grMCFMsta->GetXaxis()->SetTitle(xTitle);
     grMCFMsta->GetYaxis()->SetTitle(yTitle);

     grMCFM8sta->SetLineColor(collMCFM8);
     grMCFM8sta->SetMarkerColor(colmMCFM8);
     grMCFM8sta->SetLineWidth(1);
     grMCFM8sta->GetXaxis()->SetTitle(xTitle);
     grMCFM8sta->GetYaxis()->SetTitle(yTitle);
     grMCFM8NPsta->SetLineColor(collMCFM8);
     grMCFM8NPsta->SetMarkerColor(colmMCFM8);
     grMCFM8NPsta->SetLineWidth(1);
     grMCFM8NPsta->GetXaxis()->SetTitle(xTitle);
     grMCFM8NPsta->GetYaxis()->SetTitle(yTitle);

     grMCFM23sta->SetLineColor(collMCFM23);
     grMCFM23sta->SetMarkerColor(colmMCFM23);
     grMCFM23sta->SetLineWidth(1);
     grMCFM23sta->GetXaxis()->SetTitle(xTitle);
     grMCFM23sta->GetYaxis()->SetTitle(yTitle);
     grMCFM23NPsta->SetLineColor(collMCFM23);
     grMCFM23NPsta->SetMarkerColor(colmMCFM23);
     grMCFM23NPsta->SetLineWidth(1);
     grMCFM23NPsta->GetXaxis()->SetTitle(xTitle);
     grMCFM23NPsta->GetYaxis()->SetTitle(yTitle);

     if (setLogY && hname=="j_bjet_Pt_binned") grMCFMsta->GetYaxis()->SetRangeUser(0.01, 500.);

     if (hname=="j_2bjet_Mbb_binned") {
       grSherpafb->Draw("PCSAME");
       grMCFMsta->Draw("PCSAME");
     }
     else {
       grMCFMsta->Draw("PCSAME");
     }

     // if (hnameaMCNLO != "Z_y")  
     //   {
     // 	 /// test ---------------{
     // 	 graMCNLOfb->SetLineColor(2);
     // 	 graMCNLOfb->SetMarkerColor(2);
     // 	 graMCNLOfb->Draw("PCSAME");
     // 	 /// end test -----------}
     //   }

     grSherpafb->SetLineWidth(1);
     grSherpafb->SetLineColor(colSherpa);
     grSherpafb->SetMarkerColor(colSherpa);
     grSherpafb->Draw("PCSAME");

     grPythiafb->SetLineWidth(1);
     grPythiafb->SetLineColor(colPythia);
     grPythiafb->SetMarkerColor(colPythia);
     grPythiafb->Draw("PCSAME");

     legend->Clear();
     legend->AddEntry(grSherpafb,"SHERPA no NP","PL");
     legend->AddEntry(grPythiafb,"PYTHIA no NP","PL");
     //if (hnameaMCNLO != "Z_y")  legend->AddEntry(graMCNLOfb, "aMC\@NLO no NP","PL");	
     legend->AddEntry(grMCFMsta, "MCFM","PE");	
     legend->Draw();

     if (saveAsEps)
       {
	 std::string fname = "fig"+hname+"_MCFM_Sherpa_Pythia_noNP.eps";
	 if (setLogY) fname = "log"+fname;
	 cPlot->SaveAs(fname.c_str());
       }
     // if (saveAsRoot)
     //   {	 
     // 	 std::string fname = "fig"+hname+"_MCFM_Sherpa_Pythia_noNP.root";
     // 	 cPlot->SaveAs(fname.c_str());
     //   }

     //return;

     //if (cPlot1) delete cPlot1;
     TCanvas *cPlot1 = new TCanvas("SherpaPythia","SherpaPythia");
     if (setLogY) cPlot1->SetLogy();
     
     hempty->Draw();


     grSherpafb->SetLineStyle(2);
     grSherpafb->SetLineColor(colSherpa);
     grSherpafb->SetMarkerColor(colSherpa);
     grSherpafb->GetXaxis()->SetTitle(xTitle);
     grSherpafb->GetYaxis()->SetTitle(yTitle);
     if (setLogY && hname=="j_bjet_Pt_binned") grSherpafb->GetYaxis()->SetRangeUser(0.01, 500.);
     grSherpafb->Draw("PC");
     grPythiafb->SetLineStyle(2);
     grPythiafb->SetLineColor(colPythia);
     grPythiafb->SetMarkerColor(colPythia);
     grPythiafb->Draw("PCSAME");
     grSherpaNPfb->SetLineStyle(1);
     grSherpaNPfb->SetLineColor(colSherpa);
     grSherpaNPfb->SetMarkerColor(colSherpa);
     grSherpaNPfb->Draw("PCSAME");
     grPythiaNPfb->SetLineStyle(1);
     grPythiaNPfb->SetLineColor(colPythia);
     grPythiaNPfb->SetMarkerColor(colPythia);
     grPythiaNPfb->Draw("PCSAME");
     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO)) {
       graMCNLOMPIfb->SetLineStyle(1);
       graMCNLOMPIfb->SetLineColor(colMCNLO);
       graMCNLOMPIfb->SetMarkerColor(colMCNLO);
       graMCNLOMPIfb->Draw("PCSAME");
       if (showMCNLOnoMPI) 
       {
	   graMCNLOfb->SetLineStyle(2);
	   graMCNLOfb->SetLineColor(colMCNLO);
	   graMCNLOfb->SetMarkerColor(colMCNLO);
	   graMCNLOfb->Draw("PCSAME");
       }
     }
     legend->Clear();
     legend->AddEntry(grSherpafb,"SHERPA no NP","PEL");
     legend->AddEntry(grPythiafb,"PYTHIA no NP","PEL");
     legend->AddEntry(grSherpaNPfb,"SHERPA","PEL");
     legend->AddEntry(grPythiaNPfb,"PYTHIA","PEL");
     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))  
     {
       if (showMCNLOnoMPI) legend->AddEntry(graMCNLOfb,   "aMC@NLO no MPI","PEL");	
       legend->AddEntry(graMCNLOMPIfb, "aMC@NLO","PEL");	
     }
     legend->Draw();
     
     if (saveAsEps)
       {
	 fname = "fig"+hname+"_aMCNLO_Sherpa_Pythia.eps";
	 if (setLogY) fname = "log"+fname;
	 cPlot1->SaveAs(fname.c_str());
       }
     // if (saveAsRoot)
     //   {
     // 	 fname = "fig"+hname+"_aMCNLO_Sherpa_Pythia.root";
     // 	 cPlot1->SaveAs(fname.c_str());
     //   }     
     //return;


     //if (cPlotNP) delete cPlotNP;
     TCanvas *cPlotNP = new TCanvas("MCFMwNP","MCFMwNP");
     if (setLogY) cPlotNP->SetLogy();
     hempty->Draw();


     grMCFMNPsta_as_pdf_sca_np->SetFillColor(colfMCFMall);
     grMCFMNPsta_as_pdf_sca_np->GetXaxis()->SetTitle(xTitle);
     grMCFMNPsta_as_pdf_sca_np->GetYaxis()->SetTitle(yTitle);
     if (setLogY && hname=="j_bjet_Pt_binned") grMCFMNPsta_as_pdf_sca_np->GetYaxis()->SetRangeUser(0.01, 500.);
     grMCFMNPsta_as_pdf_sca_np->Draw("PE2");
     grMCFMNPsta_as_pdf_sca->SetFillColor(colfMCFMas_pdf_sca);
     grMCFMNPsta_as_pdf_sca->Draw("PE2SAME");
     grMCFMNPsta_as_pdf->SetFillColor(colfMCFMas_pdf);
     grMCFMNPsta_as_pdf->Draw("PE2SAME");
     grMCFMNPsta_as->SetFillColor(colfMCFMas);
     grMCFMNPsta_as->Draw("PE2SAME");
     grMCFMNPsta->SetLineColor(collMCFM);
     grMCFMNPsta->SetMarkerColor(colmMCFM);
     grMCFMNPsta->Draw("PESAME");
     // central value with mstw & NP+QED {
     grMCFM8NPsta ->Draw("PECSAME");
     grMCFM23NPsta->Draw("PECSAME");
     // central value with mstw & NP+QED }

     grSherpaNPfb->SetLineColor(colSherpa);
     grSherpaNPfb->SetLineWidth(1);
     grSherpaNPfb->SetMarkerColor(colSherpa);
     grSherpaNPfb->Draw("PECSAME");
     grPythiaNPfb->SetLineWidth(1);
     grPythiaNPfb->SetLineColor(colPythia);
     grPythiaNPfb->SetMarkerColor(colPythia);
     grPythiaNPfb->Draw("PECSAME");
     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO)) {
	 //       graMCNLOMPIfb->SetMarkerSize();
       graMCNLOMPIfb->SetLineWidth(1);
       graMCNLOMPIfb->SetLineColor(colMCNLO);
       graMCNLOMPIfb->SetMarkerColor(colMCNLO);
       graMCNLOMPIfb->Draw("PECSAME");
       // graMCNLOfb->SetLineWidth(1);
       // graMCNLOfb->SetLineStyle(2);
       // graMCNLOfb->SetLineColor(colMCNLO);
       // graMCNLOfb->SetMarkerColor(colMCNLO);
       // graMCNLOfb->Draw("PECSAME");

     }
     legend1->Clear();
     legend1->AddEntry(grMCFMNPsta_as_pdf_sca_np,"MCFM (ct10) x NP    np err.","F");
     legend1->AddEntry(grMCFMNPsta_as_pdf_sca,   "MCFM (ct10) x NP scale err.","F");
     legend1->AddEntry(grMCFMNPsta_as_pdf    ,   "MCFM (ct10) x NP   pdf err.","F");
     legend1->AddEntry(grMCFMNPsta_as        ,   "MCFM (ct10) x NP #alpha_{s} err.","F");
     legend1->AddEntry(grMCFMNPsta           ,   "MCFM (ct10) x NP stat. err.","PEL");
     legend1->AddEntry(grMCFM8NPsta          ,   "MCFM (mstw2008) x NP stat. err.","PEL");
     legend1->AddEntry(grMCFM23NPsta         ,   "MCFM (nnpdf23)  x NP stat. err.","PEL");
     legend1->AddEntry(grSherpaNPfb,"SHERPA","PEL");
     legend1->AddEntry(grPythiaNPfb,"PYTHIA","PEL");
     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))  legend1->AddEntry(graMCNLOMPIfb, "aMC@NLO","PEL");	
     //if (hnameaMCNLO != "Z_y")  legend1->AddEntry(graMCNLOfb,    "aMC@NLO no MPI","PEL");	
     legend1->Draw();
     


     if (saveAsEps)
       {
	 fname = "fig"+hname+"_MCFM_aMCNLO_Sherpa_Pythia.eps";
	 if (setLogY) fname = "log"+fname;
	 cPlotNP->SaveAs(fname.c_str());
       }
     
     //if (cPlotmcfm) delete cPlotmcfm;
     TCanvas *cPlotmcfm = new TCanvas("MCFM2pdf","MCFM2pdf");
     if (setLogY) cPlotmcfm->SetLogy();
     hempty->Draw();


     grMCFM8NPsta_as_pdf_sca_np->SetFillColor(colfMCFMall);
     grMCFM8NPsta_as_pdf_sca_np->GetXaxis()->SetTitle(xTitle);
     grMCFM8NPsta_as_pdf_sca_np->GetYaxis()->SetTitle(yTitle);
     if (setLogY && hname=="j_bjet_Pt_binned") grMCFM8NPsta_as_pdf_sca_np->GetYaxis()->SetRangeUser(0.01, 500.);
     grMCFM8NPsta_as_pdf_sca_np->Draw("PE2");
     grMCFM8NPsta_as_pdf_sca->SetFillColor(colfMCFMas_pdf_sca);
     grMCFM8NPsta_as_pdf_sca->Draw("PE2SAME");
     grMCFM8NPsta_as_pdf->SetFillColor(colfMCFMas_pdf);
     grMCFM8NPsta_as_pdf->Draw("PE2SAME");
     grMCFM8NPsta_as->SetFillColor(colfMCFMas);
     grMCFM8NPsta_as->Draw("PE2SAME");
     grMCFM8NPsta->SetLineColor(collMCFM8);
     grMCFM8NPsta->SetMarkerColor(collMCFM8);
     grMCFM8NPsta->Draw("PESAME");

     // CT10 central value with NP 
     grMCFMNPsta->Draw("PESAME");
     // MSTW2008 w and w/o NP
//      grMCFM8sta->SetLineStyle(1);
//      grMCFMsta-> SetLineStyle(2);
//      grMCFM8sta->Draw("CSAME");
//      grMCFMsta-> Draw("CSAME");

     legend2->Clear();
     legend2->AddEntry(grMCFM8NPsta_as_pdf_sca_np,"MCFM\_mstw2008 x NP    np err.","F");
     legend2->AddEntry(grMCFM8NPsta_as_pdf_sca,   "MCFM\_mstw2008 x NP scale err.","F");
     legend2->AddEntry(grMCFM8NPsta_as_pdf    ,   "MCFM\_mstw2008 x NP   pdf err.","F");
     legend2->AddEntry(grMCFM8NPsta_as        ,   "MCFM\_mstw2008 x NP #alpha_{s} err.","F");
     legend2->AddEntry(grMCFM8NPsta           ,   "MCFM\_mstw2008 x NP stat. err.","PEL");
     legend2->AddEntry(grMCFMNPsta            ,   "MCFM\_ct10 x NP stat. err.","PEL");
//      legend2->AddEntry(grMCFM8sta             ,   "MCFM\_mstw2008 stat. err.","L");
//      legend2->AddEntry(grMCFMsta              ,   "MCFM\_ct10     stat. err.","L");
     legend2->Draw();
     grMCFMsta-> SetLineStyle(1);
     

     if (saveAsRoot)
       {
	      grMCFMNPsta_as_pdf_sca_np->Write("grMCFMNPsta_as_pdf_sca_np");
	      grMCFMNPsta_as_pdf_sca   ->Write("grMCFMNPsta_as_pdf_sca");
	      grMCFMNPsta_as_pdf       ->Write("grMCFMNPsta_as_pdf");
	      grMCFMNPsta_as           ->Write("grMCFMNPsta_as");
	      grMCFMNPsta              ->Write("grMCFMNPsta");
	      grMCFM8NPsta_as_pdf_sca_np->Write("grMCFM8NPsta_as_pdf_sca_np");
	      grMCFM8NPsta_as_pdf_sca   ->Write("grMCFM8NPsta_as_pdf_sca");
	      grMCFM8NPsta_as_pdf       ->Write("grMCFM8NPsta_as_pdf");
	      grMCFM8NPsta_as           ->Write("grMCFM8NPsta_as");
	      grMCFM8NPsta              ->Write("grMCFM8NPsta");
	      grMCFM23NPsta_as_pdf_sca_np->Write("grMCFM23NPsta_as_pdf_sca_np");
	      grMCFM23NPsta_as_pdf_sca   ->Write("grMCFM23NPsta_as_pdf_sca");
	      grMCFM23NPsta_as_pdf       ->Write("grMCFM23NPsta_as_pdf");
	      grMCFM23NPsta_as           ->Write("grMCFM23NPsta_as");
	      grMCFM23NPsta              ->Write("grMCFM23NPsta");
	 //	 fname = "fig"+hname+"_MCFM_aMCNLO_Sherpa_Pythia.root";
	 //	 cPlotNP->SaveAs(fname.c_str());
       }

     if (saveAsEps)
       {
	 fname = "fig"+hname+"_MCFM2pdf.eps";
	 if (setLogY) fname = "log"+fname;
	 cPlotmcfm->SaveAs(fname.c_str());
       }
     
     //if (cPlotmcfm23) delete cPlotmcfm23;
     TCanvas *cPlotmcfm23 = new TCanvas("MCFM23pdf","MCFM23pdf");
     if (setLogY) cPlotmcfm23->SetLogy();
     hempty->Draw();


     grMCFM23NPsta_as_pdf_sca_np->SetFillColor(colfMCFMall);
     grMCFM23NPsta_as_pdf_sca_np->GetXaxis()->SetTitle(xTitle);
     grMCFM23NPsta_as_pdf_sca_np->GetYaxis()->SetTitle(yTitle);
     if (setLogY && hname=="j_bjet_Pt_binned") grMCFM23NPsta_as_pdf_sca_np->GetYaxis()->SetRangeUser(0.01, 500.);
     grMCFM23NPsta_as_pdf_sca_np->Draw("PE2");
     grMCFM23NPsta_as_pdf_sca->SetFillColor(colfMCFMas_pdf_sca);
     grMCFM23NPsta_as_pdf_sca->Draw("PE2SAME");
     grMCFM23NPsta_as_pdf->SetFillColor(colfMCFMas_pdf);
     grMCFM23NPsta_as_pdf->Draw("PE2SAME");
     grMCFM23NPsta_as->SetFillColor(colfMCFMas);
     grMCFM23NPsta_as->Draw("PE2SAME");
     //     grMCFM23NPsta->SetLineColor(collMCFM23);
     //     grMCFM23NPsta->SetMarkerColor(collMCFM23);
     grMCFM23NPsta->Draw("PESAME");

     // CT10 central value with NP 
     grMCFMNPsta->Draw("PESAME");
     // CT10 central value with NP 
     grMCFM8NPsta->Draw("PESAME");
     // MSTW2008 w and w/o NP
     /// here without NP corrections {
//      grMCFM23sta->SetLineStyle(1);
//      grMCFM8sta->SetLineStyle(3);
//      grMCFMsta-> SetLineStyle(2);
//      grMCFM23sta->Draw("CSAME");
//      grMCFM8sta-> Draw("CSAME");
//      grMCFMsta->  Draw("CSAME");
     ///  here without NP corrections }


     legend2->Clear();
     legend2->AddEntry(grMCFM23NPsta_as_pdf_sca_np,"MCFM\_nnpdf23 x NP    np err.","F");
     legend2->AddEntry(grMCFM23NPsta_as_pdf_sca,   "MCFM\_nnpdf23 x NP scale err.","F");
     legend2->AddEntry(grMCFM23NPsta_as_pdf    ,   "MCFM\_nnpdf23 x NP   pdf err.","F");
     legend2->AddEntry(grMCFM23NPsta_as        ,   "MCFM\_nnpdf23 x NP #alpha_{s} err.","F");
     legend2->AddEntry(grMCFM23NPsta           ,   "MCFM\_nnpdf23 x NP stat. err.","PEL");
     legend2->AddEntry(grMCFMNPsta            ,   "MCFM\_ct10 x NP stat. err.","PEL");
     legend2->AddEntry(grMCFM8NPsta           ,   "MCFM\_mstw2008 x NP stat. err.","PEL");
//      legend2->AddEntry(grMCFM23sta            ,   "MCFM\_nnpdf23  stat. err.","L");
//      legend2->AddEntry(grMCFMsta              ,   "MCFM\_ct10     stat. err.","L");
//      legend2->AddEntry(grMCFM8sta             ,   "MCFM\_mstw2008 stat. err.","L");
     legend2->Draw();
     grMCFMsta-> SetLineStyle(1);
     

     if (saveAsEps)
       {
	 fname = "fig"+hname+"_MCFM23pdf.eps";
	 if (setLogY) fname = "log"+fname;
	 cPlotmcfm23->SaveAs(fname.c_str());
       }
     // if (saveAsRoot)
     //   {
     // 	 fname = "fig"+hname+"_MCFM_aMCNLO_Sherpa_Pythia.root";
     // 	 cPlotNP->SaveAs(fname.c_str());
     //   }


     //if (cPlotNP_01) delete cPlotNP_01;
     TCanvas *cPlotNP_01 = new TCanvas("MCFM_w_wo_NP","MCFM_w_wo_NP");
     if (setLogY) cPlotNP_01->SetLogy();
     hempty->Draw();

//      grMCFMsta->SetLineColor(1);
//      grMCFMsta->GetXaxis()->SetTitle(xTitle);
//      grMCFMsta->GetYaxis()->SetTitle(yTitle);
//      if (setLogY && hname=="j_bjet_Pt_binned") grMCFMsta->GetYaxis()->SetRangeUser(0.01, 500.);
//      grMCFMsta->Draw("C");
     grMCFMNPsta_as_pdf_sca_np->SetFillColor(5);
     //grMCFMNPsta_as_pdf_sca_np->Draw("E2SAME");
     grMCFMNPsta_as_pdf_sca_np->Draw("E2");
     //grMCFMsta->Draw("CSAME");
     grMCFM8NPsta_as_pdf_sca_np->SetFillStyle(3004);
     grMCFM8NPsta_as_pdf_sca_np->SetFillColor(collMCFM8);
     grMCFM8NPsta_as_pdf_sca_np->Draw("E2SAME");
     //grMCFM8sta->Draw("CSAME");
     grMCFM23NPsta_as_pdf_sca_np->SetFillStyle(3005);
     grMCFM23NPsta_as_pdf_sca_np->SetFillColor(collMCFM23);
     grMCFM23NPsta_as_pdf_sca_np->Draw("E2SAME");
     //grMCFM23sta->Draw("CSAME");

     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))  
       {     
	 /// test ---------------{
	 graMCNLOMPIfb->SetLineColor(2);
	 graMCNLOMPIfb->SetMarkerColor(2);
	 graMCNLOMPIfb->Draw("PSAME");
	 /// end test -----------}
	 /// test ---------------{
	 if (showMCNLOnoMPI) 
	 {
	     graMCNLOfb->SetLineColor(2);	
	     graMCNLOfb->SetLineStyle(2);
	     graMCNLOfb->SetMarkerColor(2);	
	     graMCNLOfb->SetMarkerStyle(24);
	     graMCNLOfb->Draw("PCSAME");
	 }
	 /// end test -----------}
       }
     legend->Clear();
     legend->AddEntry(grMCFMNPsta_as_pdf_sca_np,   "MCFM (ct10) x NP",    "F");
     legend->AddEntry(grMCFM8NPsta_as_pdf_sca_np,  "MCFM (mstw2008) x NP","F");
     legend->AddEntry(grMCFM23NPsta_as_pdf_sca_np, "MCFM (nnpdf23) x NP", "F");
     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))   
       {
	 legend->AddEntry(graMCNLOMPIfb,            "aMC@NLO","PEL");
	 if (showMCNLOnoMPI) legend->AddEntry(graMCNLOfb,               "aMC@NLO no MPI","PEL");
       }
     //     legend->AddEntry(grMCFMsta           ,     "MCFM (ct10) stat. err.","L");
     //     legend->AddEntry(grMCFM8sta           ,     "MCFM (mstw2008) stat. err.","L");
     //     legend->AddEntry(grMCFM23sta           ,     "MCFM (nnpdf23) stat. err.","L");
     legend->Draw();
     


     if (saveAsEps)
       {
	 fname = "fig"+hname+"_MCFM_w_wo_NP.eps";
	 if (setLogY) fname = "log"+fname;
	 cPlotNP_01->SaveAs(fname.c_str());
       }
     // if (saveAsRoot)
     //   {
     // 	 fname = "fig"+hname+"_MCFM_w_wo_NP.root";
     // 	 cPlotNP_01->SaveAs(fname.c_str());
     //   }
     
     if (hname.substr(0,6)=="j_bjet")
       {
     
	 //if (cPlotNP_01data) delete cPlotNP_01data;
	 TCanvas *cPlotNP_01data = new TCanvas("MCFM_w_wo_NP_andData","MCFM_w_wo_NP_andData");
	 if (setLogY) cPlotNP_01data->SetLogy();
	 hempty->Draw();

	 //      grMCFMsta->SetLineColor(1);
	 //      grMCFMsta->GetXaxis()->SetTitle(xTitle);
	 //      grMCFMsta->GetYaxis()->SetTitle(yTitle);
	 //      if (setLogY && hname=="j_bjet_Pt_binned") grMCFMsta->GetYaxis()->SetRangeUser(0.01, 500.);
	 //      grMCFMsta->Draw("C");
	 grMCFMNPsta_as_pdf_sca_np->SetFillColor(5);
	 //grMCFMNPsta_as_pdf_sca_np->Draw("E2SAME");
	 grMCFMNPsta_as_pdf_sca_np->Draw("E2");
	 //grMCFMsta->Draw("CSAME");
	 grMCFM8NPsta_as_pdf_sca_np->SetFillStyle(3004);
	 grMCFM8NPsta_as_pdf_sca_np->SetFillColor(collMCFM8);
	 grMCFM8NPsta_as_pdf_sca_np->Draw("E2SAME");
	 //grMCFM8sta->Draw("CSAME");
	 grMCFM23NPsta_as_pdf_sca_np->SetFillStyle(3005);
	 grMCFM23NPsta_as_pdf_sca_np->SetFillColor(collMCFM23);
	 grMCFM23NPsta_as_pdf_sca_np->Draw("E2SAME");
	 //grMCFM23sta->Draw("CSAME");

	 if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))   
	   {     
	     /// test ---------------{
	     graMCNLOMPIfb->SetLineColor(2);
	     graMCNLOMPIfb->SetMarkerColor(2);
	     graMCNLOMPIfb->Draw("PSAME");
	     /// end test -----------}
	     /// test ---------------{
	     if (showMCNLOnoMPI)
	       {
		 graMCNLOfb->SetLineColor(2);
		 graMCNLOfb->SetLineStyle(2);
		 graMCNLOfb->SetMarkerColor(2);
		 graMCNLOfb->SetMarkerStyle(24);
		 graMCNLOfb->Draw("PCSAME");
	       }
	     /// end test -----------}
	   }

	 grdataNbbb->SetMarkerColor(1);
	 grdataNBayes->SetMarkerColor(1);
	 grdataMbbb->SetMarkerColor(1);
	 grdataMBayes->SetMarkerColor(1);
	 grdataNbbb->SetMarkerStyle(24);
	 grdataNbbb->SetMarkerSize(1.5);
	 grdataNBayes->SetMarkerStyle(20);
	 grdataMbbb->SetMarkerStyle(26);
	 grdataMbbb->SetMarkerSize(1.5);
	 grdataMBayes->SetMarkerStyle(23);
    

	 grdataMbbb->Draw("PESAME");
	 grdataMBayes->Draw("PESAME");     
	 if (hname=="j_bjet_Pt_binned"){
	   grdataNbbb->Draw("PESAME");
	   grdataNBayes->Draw("PESAME");
	 }
	 legend1->Clear();
	 legend1->AddEntry(grMCFMNPsta_as_pdf_sca_np,"MCFM (ct10) x NP","F");
	 //     legend->AddEntry(grMCFMsta           ,     "MCFM (ct10) stat. err.","L");
	 legend1->AddEntry(grMCFM8NPsta_as_pdf_sca_np,"MCFM (mstw2008) x NP","F");
	 //     legend->AddEntry(grMCFM8sta           ,     "MCFM (mstw2008) stat. err.","L");
	 legend1->AddEntry(grMCFM23NPsta_as_pdf_sca_np,"MCFM (nnpdf23) x NP","F");
	 //     legend->AddEntry(grMCFM23sta           ,     "MCFM (nnpdf23) stat. err.","L");
	 if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))  
	   {
	     legend1->AddEntry(graMCNLOMPIfb,            "aMC@NLO","PEL");
	     if (showMCNLOnoMPI) legend1->AddEntry(graMCNLOfb,               "aMC@NLO no MPI","PEL");
	     legend1->AddEntry(grdataMbbb,               "Mark bxb Alpgen-unf.","PEL");
	     legend1->AddEntry(grdataMBayes,             "Mark Bayes Alpgen-unf.","PEL");
	     if (hname=="j_bjet_Pt_binned"){
	       legend1->AddEntry(grdataNbbb,               "Nicola bxb Alpgen-unf.","PEL");
	       legend1->AddEntry(grdataNBayes,             "Nicola Bayes Alpgen-unf.","PEL");
	     }
	   }
	 legend1->Draw();
     


	 if (saveAsEps)
	   {
	     fname = "fig"+hname+"_MCFM_w_wo_NP_data.eps";
	     if (setLogY) fname = "log"+fname;
	     cPlotNP_01data->SaveAs(fname.c_str());
	   }
	 // if (saveAsRoot)
	 //   {
	 // 	 fname = "fig"+hname+"_MCFM_w_wo_NP.root";
	 // 	 cPlotNP_01->SaveAs(fname.c_str());
	 //   }

     }
     
     if (saveAsRoot) 
       {
	 foutRoot->Close();
       }



     // TCanvas *cPlotNP1 = new TCanvas("cPlotNP1","myDistributionNP1");
     // grSherpa->SetLineColor(3);
     // grSherpa->SetMarkerColor(3);
     // grSherpa->Draw("AP");
     // grPythia->SetLineColor(4);
     // grPythia->SetMarkerColor(4);
     // grPythia->Draw("PSAME");
     // graMCNLOMPI->SetLineColor(2);
     // graMCNLOMPI->SetMarkerColor(2);
     // graMCNLOMPI->Draw("PSAME");


}
/*
 TFile**		../MCFM/Zb/ZbMCFMgraphCT10.root	
 TFile*		../MCFM/Zb/ZbMCFMgraphCT10.root	
  KEY: TGraphAsymmErrors	Pt_gae_stat;1	Graph
  KEY: TGraphAsymmErrors	Pt_gae_scale;1	Graph
  KEY: TGraphAsymmErrors	Pt_gae_pdf;1	Graph
  KEY: TGraphAsymmErrors	Pt_gae_as;1	Graph
  KEY: TGraphAsymmErrors	Y_gae_stat;1	Graph
  KEY: TGraphAsymmErrors	Y_gae_scale;1	Graph
  KEY: TGraphAsymmErrors	Y_gae_pdf;1	Graph
  KEY: TGraphAsymmErrors	Y_gae_as;1	Graph
  KEY: TGraphAsymmErrors	ZPt_gae_stat;1	Graph
  KEY: TGraphAsymmErrors	ZPt_gae_scale;1	Graph
  KEY: TGraphAsymmErrors	ZPt_gae_pdf;1	Graph
  KEY: TGraphAsymmErrors	ZPt_gae_as;1	Graph
  KEY: TGraphAsymmErrors	ZY_gae_stat;1	Graph
  KEY: TGraphAsymmErrors	ZY_gae_scale;1	Graph
  KEY: TGraphAsymmErrors	ZY_gae_pdf;1	Graph
  KEY: TGraphAsymmErrors	ZY_gae_as;1	Graph
  KEY: TGraphAsymmErrors	ZdR_gae_stat;1	Graph
  KEY: TGraphAsymmErrors	ZdR_gae_scale;1	Graph
  KEY: TGraphAsymmErrors	ZdR_gae_pdf;1	Graph
  KEY: TGraphAsymmErrors	ZdR_gae_as;1	Graph
  KEY: TGraphAsymmErrors	ZdPhi_gae_stat;1	Graph
  KEY: TGraphAsymmErrors	ZdPhi_gae_scale;1	Graph
  KEY: TGraphAsymmErrors	ZdPhi_gae_pdf;1	Graph
  KEY: TGraphAsymmErrors	ZdPhi_gae_as;1	Graph
*/
/*
 TFile*		../aMCatNLO/merged.root	 // contains NP 
  OBJ: TH1F	Z_pT	 : 0 at: 0x18c0200
  OBJ: TH1F	Z_pT_bb	 : 0 at: 0x1bdcc00
  KEY: TH1F	Z_pT;1	
  KEY: TH1F	Zj_dphi;1	
  KEY: TH1F	Zj_dR;1	
  KEY: TH1F	bjet_pT;1	
  KEY: TH1F	bjet_y;1	
  KEY: TH1F	Z_pT_bb;1	
  KEY: TH1F	bb_mass;1	
  KEY: TH1F	bb_dR;1	
TFile**		../MCFM/Zbb/ZbbMCFMgraphCT10.root	
 TFile*		../MCFM/Zbb/ZbbMCFMgraphCT10.root	
  KEY: TGraphAsymmErrors	ZPT_gae_stat;1	Graph
  KEY: TGraphAsymmErrors	ZPT_gae_scale;1	Graph
  KEY: TGraphAsymmErrors	ZPT_gae_pdf;1	Graph
  KEY: TGraphAsymmErrors	ZPT_gae_as;1	Graph
  KEY: TGraphAsymmErrors	bbdR_gae_stat;1	Graph
  KEY: TGraphAsymmErrors	bbdR_gae_scale;1	Graph
  KEY: TGraphAsymmErrors	bbdR_gae_pdf;1	Graph
  KEY: TGraphAsymmErrors	bbdR_gae_as;1	Graph
  KEY: TGraphAsymmErrors	bbM_gae_stat;1	Graph
  KEY: TGraphAsymmErrors	bbM_gae_scale;1	Graph
  KEY: TGraphAsymmErrors	bbM_gae_pdf;1	Graph
  KEY: TGraphAsymmErrors	bbM_gae_as;1	Graph
*/


/*****.....**** Nicola yield 
17889.2
16641.7
7705.4
3704.6
1720.5
300.4
*///////////////
/***********  Nicola rooUnfold
17694.4
16623
7726.9
3736.4
1705.7
301.8
*///////////////
/***** Mark bin by bin Alpgen 
21230.8
16282.5
7007.6
3292
1493.2
259.3
*///////////////
/***** Mark bin by bin Sherpa
21870.3
16553
7299.7
3488.9
1675.4
294.1
*///////////////