#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"
#include "getGraphRatio.C"

#define saveAsRoot      false
#define saveAsEps       true
#define showMCNLOnoMPI  false
#define m_debug         0
#define writeDataResults true
#define MBayes true
#define NBayes false
#define Mbbb   false
#define Nbbb   false

// TGraphAsymmErrors* getGraphRatio(TGraphAsymmErrors* grD, TGraphAsymmErrors* grN)
// {
//   Double_t* grDY = grD->GetY();
//   Double_t* grNY = grN->GetY();
//   Double_t* grNEYl = grN->GetEYlow();
//   Double_t* grNEYh = grN->GetEYhigh();
//   Double_t ratY[100];
//   Double_t ratEYl[100];
//   Double_t ratEYh[100];

//   for (int k=0; k<grD->GetN(); ++k)
//     {
//       ratY[k] =0.;
//       ratEYl[k] =0.;
//       ratEYh[k] =0.;
//       if (grDY[k]!=0.) {
// 	ratY[k] = grNY[k]/grDY[k];
// 	ratEYl[k] = grNEYl[k]/grDY[k];
// 	ratEYh[k] = grNEYh[k]/grDY[k];
//       }
//     }

//   TGraphAsymmErrors* grR = new TGraphAsymmErrors(grD->GetN(), grD->GetX(), ratY, grD->GetEXlow(), grD->GetEXhigh(), ratEYl, ratEYh);
//   grR->SetLineStyle(grN->GetLineStyle());
//   grR->SetLineColor(grN->GetLineColor());
//   grR->SetLineWidth(grN->GetLineWidth());
//   //  grR->SetLineAttribute(grN->GetLineAttribute());

//   grR->SetMarkerStyle(grN->GetMarkerStyle());
//   grR->SetMarkerColor(grN->GetMarkerColor());
//   //  grR->SetMarkerWidth(grN->GetMarkerWidth());
//   //  grR->SetMarkerAttribute(grN->GetMarkerAttribute());

//   grR->SetFillStyle(grN->GetFillStyle());
//   grR->SetFillColor(grN->GetFillColor());
//   //  grR->SetFillWidth(grN->GetFillWidth());
//   //  grR->SetFillAttribute(grN->GetFillAttribute());

//   return grR;
// }

void makeDistribPlotWratio(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; ////// value used  for v. 1.4
  ////////  Double_t Lumi = 4.727061 // fb  --- according to 4.727
  // Lumi = 4.57944; ///// new - latest round 
  Double_t kfactor=1.25;
  

  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* fileaMCNLO4f = new TFile("/afs/le.infn.it/user/s/spagnolo/atlas/Physics/Generators/aMCatNLO/Chiodini/resultsWithRivet/Rivet_4fHerwPP_10E4.root");
//      TFile* fileaMCNLO4f = new TFile("/afs/le.infn.it/project/atlasewk/Generators/aMC_at_NLO/rootFilesUCL/cnvHepMC/Zbb_all.v3.Rivet.root");
//      TFile* fileaMCNLO4f_v3 = new TFile("/afs/le.infn.it/project/atlasewk/Generators/aMC_at_NLO/lecceProd/oldsetup13Aug/histo_oldSet4fZbb.v3.Rivet.root");
//      TFile* fileaMCNLO5fZb_v3 = new TFile("/afs/le.infn.it/project/atlasewk/Generators/aMC_at_NLO/lecceProd/oldsetup13Aug/histo_oldSet5fZb.v3.Rivet.root");
//      TFile* fileaMCNLO5fZb_MI_v3 = new TFile("/afs/le.infn.it/project/atlasewk/Generators/aMC_at_NLO/lecceProd/oldsetup13Aug/histo_oldSet5fZb_MI.v3.Rivet.root");
//      TFile* fileaMCNLO5fZb_noFR_v3 = new TFile("/afs/le.infn.it/project/atlasewk/Generators/aMC_at_NLO/lecceProd/oldsetup13Aug/histo_oldSet5fZb_noFR.v3.Rivet.root");
//      TFile* fileaMCNLO5fZj_v3 = new TFile("/afs/le.infn.it/project/atlasewk/Generators/aMC_at_NLO/lecceProd/oldsetup13Aug/histo_oldSet5fZj.v3.Rivet.root");
     TFile* fileaMCNLO4f           = new TFile("../aMCatNLO/Zbb_all.v3.Rivet.root");
     TFile* fileaMCNLO4f_v3        = new TFile("../aMCatNLO/histo_oldSet4fZbb.v3.Rivet.root");
     TFile* fileaMCNLO5fZb_v3      = new TFile("../aMCatNLO/histo_oldSet5fZb.v3.Rivet.root");
     TFile* fileaMCNLO5fZb_MI_v3   = new TFile("../aMCatNLO/histo_oldSet5fZb_MI.v3.Rivet.root");
     TFile* fileaMCNLO5fZb_noFR_v3 = new TFile("../aMCatNLO/histo_oldSet5fZb_noFR.v3.Rivet.root");
     TFile* fileaMCNLO5fZj_v3      = new TFile("../aMCatNLO/histo_oldSet5fZj.v3.Rivet.root");
     //     TFile* fileaMCNLO4f = new TFile("/afs/le.infn.it/user/s/spagnolo/atlas/Physics/Generators/aMCatNLO/Chiodini/resultsWithRivet/Rivet_4fHerwPP_10E4.root");
     TFile* fileMCFM;
     TFile* fileMCFMmstw;
     TFile* fileMCFMnnpdf;
     TFile* fileAlpgen;
     if (hname.substr(0,3)=="j_b") 
     {
       fileMCFM      = new TFile("../MCFM/Zb/ZbMCFMgraphCT10excl266.root");
       fileMCFMmstw  = new TFile("../MCFMmstw/Zb/ZbMCFMgraphMSTW2008excl266.root");
       fileMCFMnnpdf = new TFile("../MCFMnnpdf/Zb/ZbMCFMgraphNNPDF23excl266.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");
       fileAlpgen = 0;
     }



     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 *graMCNLO4f; 
     TGraphAsymmErrors *graMCNLO4fls; 
     TGraphAsymmErrors *graMCNLO5fZb; 
     TGraphAsymmErrors *graMCNLO5fZj; 
     TGraphAsymmErrors *graMCNLO5fZb_MI; 
     TGraphAsymmErrors *graMCNLO5fZb_noFR; 
     TGraphAsymmErrors *grSherpafb; // in fb 
     TGraphAsymmErrors *grPythiafb; // in fb 
     TGraphAsymmErrors *graMCNLOfb; 
     TGraphAsymmErrors *graMCNLO4ffb; 
     TGraphAsymmErrors *graMCNLO4flsfb; 
     TGraphAsymmErrors *graMCNLO5fZbfb; 
     TGraphAsymmErrors *graMCNLO5fZjfb; 
     TGraphAsymmErrors *graMCNLO5fZb_MIfb; 
     TGraphAsymmErrors *graMCNLO5fZb_noFRfb; 
     // sherpa and pythia with NP 
     TGraphAsymmErrors *grAlpgenNP; 
     TGraphAsymmErrors *grSherpaNP; 
     TGraphAsymmErrors *grPythiaNP; 
     TGraphAsymmErrors *graMCNLOMPI; 
     TGraphAsymmErrors *grAlpgenNPfb; // in fb  
     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;  
     TLegend *legendDMC = 0;  
     TLegend *legendRATD = 0;  
     TLegend *legendRATS = 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;
     
     int colfMCFMstat  = 0;
     int colfMCFMas    = 5;
     int colfMCFMas_pdf  = 2;
     int colfMCFMas_pdf_sca  = 94;
     int colfMCFMall = 6;
     int colSherpa   = 3;
     int colAlpgen   = 6;
     int colPythia   = 4;
     int collMCFM    = 216;
     int colmMCFM    = 216;
     int collMCFM8   = 1;
     int colmMCFM8   = 1;
     int collMCFM23  = 7;
     int colmMCFM23  = 7;
     
     std::string fdataname = "../data/Unfolded_"+hname+"_SYST.txt";
     if (m_debug > -1) 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);
     // TGraph* grEmptyRat = new TGraph(n0,xEmpty,yEmpty);
     

     TH2F* hempty; 
     TH2F* hemptyRat; 
     TH2F* hemptyRatF; 

     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];
     Double_t edataMBayesSup[100];
     Double_t edataMBayesSdw[100];
     Double_t edataNBayesSup[100];
     Double_t edataNBayesSdw[100];
     Double_t edataMBayesSSup[100];
     Double_t edataMBayesSSdw[100];
     Double_t edataNBayesSSup[100];
     Double_t edataNBayesSSdw[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.;	 
	 edataMBayesSup[j] = 0.;	 
	 edataMBayesSdw[j] = 0.;	 
	 edataNBayesSup[j] = 0.;	 
	 edataNBayesSdw[j] = 0.;	 
	 edataMBayesSSup[j] = 0.;	 
	 edataMBayesSSdw[j] = 0.;	 
	 edataNBayesSSup[j] = 0.;	 
	 edataNBayesSSdw[j] = 0.;	 
     }
     

     double ratfmin=0.8;
     double ratfmax=1.2;
     
     if (hname=="j_bjet_Pt_binned")
       {
	 fileAlpgen= new TFile("../alpgen/UnfoldingSummary_CombinedChannel_JetPt_Central.root");
	 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
	 if (setLogY) hemptyRat = new TH2F("emptyRat","emptyRat",10,0.0,600.,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hemptyRat =         new TH2F("emptyRat","emptyRat",10,0.0,600.,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hemptyRatF =         new TH2F("emptyRat","emptyRat",10,0.0,600.,5,ratfmin,ratfmax); // 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)
       {
	 fileAlpgen= new TFile("../alpgen/UnfoldingSummary_CombinedChannel_JetY_Central.root");
	 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
	 if (setLogY) hemptyRat = new TH2F("emptyRat","emptyRat",10,0.,3.,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hemptyRat =         new TH2F("emptyRat","emptyRat",10,0.,3.,5,.4,1.6); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hemptyRatF =             new TH2F("emptyRat","emptyRat",10,0.,3.,5,ratfmin,ratfmax); // 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)
       {
	 fileAlpgen= new TFile("../alpgen/UnfoldingSummary_CombinedChannel_ZY_Central.root");
	 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
	 if (setLogY) hemptyRat = new TH2F("emptyRat","emptyRat",10,0.,3.,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hemptyRat =         new TH2F("emptyRat","emptyRat",10,0.,3.,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hemptyRatF =             new TH2F("emptyRat","emptyRat",10,0.,3.,5,ratfmin,ratfmax); // 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)
       {
	 fileAlpgen= new TFile("../alpgen/UnfoldingSummary_CombinedChannel_ZPt_Central.root");
	 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
	 if (setLogY) hemptyRat = new TH2F("emptyRat","emptyRat",10,0.,600.,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hemptyRat =         new TH2F("emptyRat","emptyRat",10,0.,600.,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hemptyRatF =             new TH2F("emptyRat","emptyRat",10,0.,600.,5,ratfmin,ratfmax); // 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)
       {
	 fileAlpgen= new TFile("../alpgen/UnfoldingSummary_CombinedChannel_DR_Central.root");
	 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,6.,10,3.,10000.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hempty = new TH2F("empty","empty",10,0.0,6.,10,0.,5500.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 if (setLogY) hemptyRat = new TH2F("emptyRat","emptyRat",10,0.,6.,5,.07,3.5); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hemptyRat =         new TH2F("emptyRat","emptyRat",10,0.,6.,5,.05,2.5); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hemptyRatF =             new TH2F("emptyRat","emptyRat",10,0.,6.,5,ratfmin,ratfmax); // 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)
       {
	 fileAlpgen= new TFile("../alpgen/UnfoldingSummary_CombinedChannel_DPhi_Central.root");
	 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
	 if (setLogY) hemptyRat = new TH2F("emptyRat","emptyRat",10,0.,3.5,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hemptyRat =         new TH2F("emptyRat","emptyRat",10,0.,3.5,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hemptyRatF =             new TH2F("emptyRat","emptyRat",10,0.,3.5,5,ratfmin,ratfmax); // 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
	 if (setLogY) hemptyRat = new TH2F("emptyRat","emptyRat",10,0.,300.,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hemptyRat =         new TH2F("emptyRat","emptyRat",10,0.,300.,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hemptyRatF =             new TH2F("emptyRat","emptyRat",10,0.,300.,5,ratfmin,ratfmax); // 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 statistic
	 if (setLogY) hemptyRat = new TH2F("emptyRat","emptyRat",10,0.,5.2,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hemptyRat =         new TH2F("emptyRat","emptyRat",10,0.,5.2,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hemptyRatF =             new TH2F("emptyRat","emptyRat",10,0.,5.2,5,ratfmin,ratfmax); // 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
	 if (setLogY) hemptyRat = new TH2F("emptyRat","emptyRat",10,0.,400.,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hemptyRat =         new TH2F("emptyRat","emptyRat",10,0.,400.,5,.3,2.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 hemptyRatF =             new TH2F("emptyRat","emptyRat",10,0.,400.,5,ratfmin,ratfmax); // 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);
     hemptyRat->GetXaxis()->SetTitle(xTitle);
     hemptyRat->GetYaxis()->SetTitle("Ratio to data");
     hemptyRatF->GetXaxis()->SetTitle(xTitle);
     hemptyRatF->GetYaxis()->SetTitle("Ratio to Analysis (1) Bayesian");
     //TH2F* hemptyCopy;


     if (legend==0) {
       if ("j_2bjet_dR_binned"==hname || "j_bjet_ZdPhiF20_binned"==hname || "j_bjet_ZdR620_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 || "j_bjet_ZdR620_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 || "j_bjet_ZdR620_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);  
     if (legendDMC==0) {
       if(  "j_bjet_ZdPhiF20_binned"==hname || "j_bjet_ZdR620_binned"==hname) legendDMC=new TLegend(0.175,0.65,0.475,0.90);
       else legendDMC=new TLegend(0.53,0.65,0.98,0.90);  
     }
     legendDMC->SetShadowColor(1);
     legendDMC->SetBorderSize(0);
     legendDMC->SetFillStyle(0);
     //legend->SetFillColor(kWhite);
     legendDMC->SetTextSize(0.04);  
     if (legendRATD==0) {
       //if( "j_bjet_ZdR620_binned"==hname) legendRATD=new TLegend(0.175,0.55,0.475,0.90);
       legendRATD=new TLegend(0.17,0.80,0.5,0.97);  
     }
     legendRATD->SetShadowColor(1);
     legendRATD->SetBorderSize(0);
     legendRATD->SetFillStyle(0);
     //legend->SetFillColor(kWhite);
     legendRATD->SetTextSize(0.04);  
     if (legendRATS==0) {
       //legendRATS=new TLegend(0.5,0.70,0.98,0.97);   // with 4 theories 
       legendRATS=new TLegend(0.5,0.80,0.98,0.97);   // with 2 theories 
     }
     legendRATS->SetShadowColor(1);
     legendRATS->SetBorderSize(0);
     legendRATS->SetFillStyle(0);
     //legend->SetFillColor(kWhite);
     legendRATS->SetTextSize(0.04);  
     
     // 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());
     if (fileAlpgen!=NULL) {
	 TH1* htmp= NULL;
	 htmp =  (TH1*)fileAlpgen->Get("h_MC");
	 //if (htmp != NULL) grAlpgenNP = new TGraphAsymmErrors(htmp);
	 grAlpgenNP = (TGraphAsymmErrors*)fileAlpgen->Get("gmc");
     }
     else grAlpgenNP=NULL;
     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());
     graMCNLO4f = (TGraphAsymmErrors*)fileaMCNLO4f->Get(hname.c_str());
     graMCNLO4fls = (TGraphAsymmErrors*)fileaMCNLO4f_v3->Get(hname.c_str());
     graMCNLO5fZb      = (TGraphAsymmErrors*)fileaMCNLO5fZb_v3     ->Get(hname.c_str());
     graMCNLO5fZj      = (TGraphAsymmErrors*)fileaMCNLO5fZj_v3     ->Get(hname.c_str());
     graMCNLO5fZb_MI   = (TGraphAsymmErrors*)fileaMCNLO5fZb_MI_v3  ->Get(hname.c_str());
     graMCNLO5fZb_noFR = (TGraphAsymmErrors*)fileaMCNLO5fZb_noFR_v3->Get(hname.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  grAlpgenNPYfb[100]; 
     Double_t  grAlpgenNPEyLowfb[100];
     Double_t  grAlpgenNPEyHighfb[100];
     Double_t  grAlpgenNPYfb_nokf[100]; 
     Double_t  grAlpgenNPEyLowfb_nokf[100];
     Double_t  grAlpgenNPEyHighfb_nokf[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* graMCNLO4fY;        graMCNLO4fY      = graMCNLO4f->GetY();
     Double_t* graMCNLO4fEyLow;    graMCNLO4fEyLow  = graMCNLO4f->GetEYlow();
     Double_t* graMCNLO4fEyHigh;   graMCNLO4fEyHigh = graMCNLO4f->GetEYhigh();
     Double_t* graMCNLO4flsY;        graMCNLO4flsY      = graMCNLO4fls->GetY();
     Double_t* graMCNLO4flsEyLow;    graMCNLO4flsEyLow  = graMCNLO4fls->GetEYlow();
     Double_t* graMCNLO4flsEyHigh;   graMCNLO4flsEyHigh = graMCNLO4fls->GetEYhigh();
     Double_t* graMCNLO5fZbY;        graMCNLO5fZbY      = graMCNLO5fZb->GetY();
     Double_t* graMCNLO5fZbEyLow;    graMCNLO5fZbEyLow  = graMCNLO5fZb->GetEYlow();
     Double_t* graMCNLO5fZbEyHigh;   graMCNLO5fZbEyHigh = graMCNLO5fZb->GetEYhigh();
     Double_t* graMCNLO5fZjY;        graMCNLO5fZjY      = graMCNLO5fZj->GetY();
     Double_t* graMCNLO5fZjEyLow;    graMCNLO5fZjEyLow  = graMCNLO5fZj->GetEYlow();
     Double_t* graMCNLO5fZjEyHigh;   graMCNLO5fZjEyHigh = graMCNLO5fZj->GetEYhigh();
     Double_t* graMCNLO5fZb_MIY;        graMCNLO5fZb_MIY      = graMCNLO5fZb_MI->GetY();
     Double_t* graMCNLO5fZb_MIEyLow;    graMCNLO5fZb_MIEyLow  = graMCNLO5fZb_MI->GetEYlow();
     Double_t* graMCNLO5fZb_MIEyHigh;   graMCNLO5fZb_MIEyHigh = graMCNLO5fZb_MI->GetEYhigh();
     Double_t* graMCNLO5fZb_noFRY;        graMCNLO5fZb_noFRY      = graMCNLO5fZb_noFR->GetY();
     Double_t* graMCNLO5fZb_noFREyLow;    graMCNLO5fZb_noFREyLow  = graMCNLO5fZb_noFR->GetEYlow();
     Double_t* graMCNLO5fZb_noFREyHigh;   graMCNLO5fZb_noFREyHigh = graMCNLO5fZb_noFR->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];
     Double_t  graMCNLO4fYfb[100]; 
     Double_t  graMCNLO4fEyLowfb[100];
     Double_t  graMCNLO4fEyHighfb[100];
     Double_t  graMCNLO4flsYfb[100]; 
     Double_t  graMCNLO4flsEyLowfb[100];
     Double_t  graMCNLO4flsEyHighfb[100];
     Double_t  graMCNLO5fZbYfb[100]; 
     Double_t  graMCNLO5fZbEyLowfb[100];
     Double_t  graMCNLO5fZbEyHighfb[100];
     Double_t  graMCNLO5fZjYfb[100]; 
     Double_t  graMCNLO5fZjEyLowfb[100];
     Double_t  graMCNLO5fZjEyHighfb[100];
     Double_t  graMCNLO5fZb_MIYfb[100]; 
     Double_t  graMCNLO5fZb_MIEyLowfb[100];
     Double_t  graMCNLO5fZb_MIEyHighfb[100];
     Double_t  graMCNLO5fZb_noFRYfb[100]; 
     Double_t  graMCNLO5fZb_noFREyLowfb[100];
     Double_t  graMCNLO5fZb_noFREyHighfb[100];

     //return;
//      // from 29 Jun 2013 
//      if (m_debug>1) std::cout<<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% reading data"<<std::endl;     
//      Double_t X1[100], X2[100], EX2[100], X3[100], X4[100];
//      Int_t nrows = 0;
//      if (hname.substr(0,6)=="j_bjet") {
//        readVectorsFromFile(fdataname,nrows,X1,X2,EX2,X3,X4,m_debug);
//        for (int j=0; j<nrows; ++j)
// 	 {
// 	   dataNbbb[j]   = X1[j];
// 	   dataNBayes[j] = X2[j];
// 	   edataNBayes[j]= EX2[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;
      // from 6 July 2013 
      if (m_debug>1) std::cout<<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% reading data"<<std::endl;     
      Double_t X1[100], X2[100], EX2[100], X3[100], X4[100], EX4[100], sUP[100], sDW[100];
      Int_t nrows = 0;
      if (hname.substr(0,6)=="j_bjet") {
	  readVectorsFromFile(fdataname,nrows,X1,X2,EX2,X3,X4,EX4,sUP,sDW,m_debug);
        for (int j=0; j<nrows; ++j)
 	 {
 	   dataNbbb[j]   = X1[j];
 	   dataNBayes[j] = X2[j];
 	   edataNBayes[j]= EX2[j];

 	   dataMbbb[j]   = X3[j];
 	   dataMBayes[j] = X4[j];
 	   edataMBayes[j] = EX4[j];

 	   edataMBayesSup[j]= sUP[j];
 	   edataMBayesSdw[j]= sDW[j];
	   
 	   edataNBayesSup[j] = 0.15*dataNBayes[j];
 	   edataNBayesSdw[j] = 0.15*dataNBayes[j];	   

 	   edataMBayesSSup[j]= sqrt( edataMBayesSup[j]*edataMBayesSup[j] + edataMBayes[j]*edataMBayes[j] );
 	   edataMBayesSSdw[j]= sqrt( edataMBayesSdw[j]*edataMBayesSdw[j] + edataMBayes[j]*edataMBayes[j] );
 	   edataNBayesSSup[j]= sqrt( edataNBayesSup[j]*edataNBayesSup[j] + edataNBayes[j]*edataNBayes[j] );
 	   edataNBayesSSdw[j]= sqrt( edataNBayesSdw[j]*edataNBayesSdw[j] + edataNBayes[j]*edataNBayes[j] );

 	   if (m_debug > 0) std::cout<<" data == stat N, Up/Dw N, stat M, Up/Dw M"<<edataNBayes[j]<<"  "<<edataNBayesSup[j]<<"/"<<edataNBayesSdw[j]<<" -------- "<<edataMBayes[j]<<"  "<<edataMBayesSup[j]<<"/"<<edataMBayesSdw[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];
	   if ( edataMBayes[0]== 0. ) edataMBayes[k]=edataNBayes[k]; // sqrt(dataMBayes[k]);
	   edataNbbb[k]=edataNBayes[k];
	   edataMbbb[k]=edataMBayes[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);

	   edataMBayesSup[k]=0.5 * edataMBayesSup[k]/(dx*Lumi);
	   edataMBayesSdw[k]=0.5 * edataMBayesSdw[k]/(dx*Lumi);
	   edataNBayesSup[k]=0.5 * edataNBayesSup[k]/(dx*Lumi);
	   edataNBayesSdw[k]=0.5 * edataNBayesSdw[k]/(dx*Lumi);

	   edataMBayesSSup[k]=0.5 * edataMBayesSSup[k]/(dx*Lumi);
	   edataMBayesSSdw[k]=0.5 * edataMBayesSSdw[k]/(dx*Lumi);
	   edataNBayesSSup[k]=0.5 * edataNBayesSSup[k]/(dx*Lumi);
	   edataNBayesSSdw[k]=0.5 * edataNBayesSSdw[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];

	 if (grAlpgenNP!=NULL) {
	   grAlpgenNPYfb[k]      = 0.5*grAlpgenNP->GetY()[k]/Lumi/dx;
	   grAlpgenNPEyLowfb[k]  = 0.5*grAlpgenNP->GetEYlow()[k]/Lumi/dx;
	   grAlpgenNPEyHighfb[k] = 0.5*grAlpgenNP->GetEYhigh()[k]/Lumi/dx;
	   grAlpgenNPYfb_nokf[k]      = 0.5*grAlpgenNP->GetY()[k]/Lumi/dx/kfactor;
	   grAlpgenNPEyLowfb_nokf[k]  = 0.5*grAlpgenNP->GetEYlow()[k]/Lumi/dx/kfactor;
	   grAlpgenNPEyHighfb_nokf[k] = 0.5*grAlpgenNP->GetEYhigh()[k]/Lumi/dx/kfactor;
	 }	  
	 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];

	 graMCNLO4fYfb[k]      = 1000.*graMCNLO4fY[k];
	 graMCNLO4fEyLowfb[k]  = 1000.*graMCNLO4fEyLow[k];
	 graMCNLO4fEyHighfb[k] = 1000.*graMCNLO4fEyHigh[k];

	 graMCNLO4flsYfb[k]      = 1000.*graMCNLO4flsY[k];
	 graMCNLO4flsEyLowfb[k]  = 1000.*graMCNLO4flsEyLow[k];
	 graMCNLO4flsEyHighfb[k] = 1000.*graMCNLO4flsEyHigh[k];

	 graMCNLO5fZbYfb[k]      = 1000.*graMCNLO5fZbY[k];
	 graMCNLO5fZbEyLowfb[k]  = 1000.*graMCNLO5fZbEyLow[k];
	 graMCNLO5fZbEyHighfb[k] = 1000.*graMCNLO5fZbEyHigh[k];

	 graMCNLO5fZjYfb[k]      = 1000.*graMCNLO5fZjY[k];
	 graMCNLO5fZjEyLowfb[k]  = 1000.*graMCNLO5fZjEyLow[k];
	 graMCNLO5fZjEyHighfb[k] = 1000.*graMCNLO5fZjEyHigh[k];

	 graMCNLO5fZb_MIYfb[k]      = 1000.*graMCNLO5fZb_MIY[k];
	 graMCNLO5fZb_MIEyLowfb[k]  = 1000.*graMCNLO5fZb_MIEyLow[k];
	 graMCNLO5fZb_MIEyHighfb[k] = 1000.*graMCNLO5fZb_MIEyHigh[k];

	 graMCNLO5fZb_noFRYfb[k]      = 1000.*graMCNLO5fZb_noFRY[k];
	 graMCNLO5fZb_noFREyLowfb[k]  = 1000.*graMCNLO5fZb_noFREyLow[k];
	 graMCNLO5fZb_noFREyHighfb[k] = 1000.*graMCNLO5fZb_noFREyHigh[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;
     TGraphAsymmErrors *grdataNBayesSyst;
     TGraphAsymmErrors *grdataNBayesTot;
     if (hname.substr(0,6)=="j_bjet") {
	 grdataNBayes     = new TGraphAsymmErrors(Npoints, 
					      grX,                             dataNBayes, 
					      grEX,                            grEX,
					      edataNBayes,                     edataNBayes);
	 grdataNBayesSyst = new TGraphAsymmErrors(Npoints, 
					      grX,                             dataNBayes, 
					      grEX,                            grEX,
					      edataNBayesSdw,                  edataNBayesSup);
	 grdataNBayesTot  = new TGraphAsymmErrors(Npoints, 
					      grX,                             dataNBayes, 
					      grEX,                            grEX,
					      edataNBayesSSdw,                 edataNBayesSSup);
     }
     
     TGraphAsymmErrors *grdataMBayes;
     TGraphAsymmErrors *grdataMBayesSyst;
     TGraphAsymmErrors *grdataMBayesTot;
     if (hname.substr(0,6)=="j_bjet") {
	 grdataMBayes     = new TGraphAsymmErrors(Npoints, 
						  grX,                             dataMBayes, 
						  grEX,                            grEX,
						  edataMBayes,                     edataMBayes);
	 grdataMBayesSyst = new TGraphAsymmErrors(Npoints, 
						  grX,                             dataMBayes, 
						  grEX,                            grEX,
						  edataMBayesSdw,                  edataMBayesSup);
	 grdataMBayesTot  = new TGraphAsymmErrors(Npoints, 
						  grX,                             dataMBayes, 
						  grEX,                            grEX,
						  edataMBayesSSdw,                 edataMBayesSSup);
     }
     
     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 *grAlpgenNPfb;
     TGraphAsymmErrors *grAlpgenNPfb_nokf;
     if (grAlpgenNP != NULL) {
	 grAlpgenNPfb = new TGraphAsymmErrors(Npoints, 
					      grX,                  grAlpgenNPYfb, 
					      grEX,                 grEX,
					      grAlpgenNPEyLowfb,    grAlpgenNPEyHighfb);
	 grAlpgenNPfb_nokf = new TGraphAsymmErrors(Npoints, 
						   grX,                       grAlpgenNPYfb_nokf, 
						   grEX,                      grEX,
						   grAlpgenNPEyLowfb_nokf,    grAlpgenNPEyHighfb_nokf);
     }     
     TGraphAsymmErrors *grPythiafb = new TGraphAsymmErrors(Npoints, 
							   grPythia->GetX(),     grPythiaYfb, 
							   grEX,                 grEX,
							   grPythiaEyLowfb,      grPythiaEyHighfb);
     TGraphAsymmErrors *grPythiaNPfb = new TGraphAsymmErrors(Npoints, 
							     grPythia->GetX(),     grPythiaNPYfb, 
							     grEX,                 grEX,
							     grPythiaNPEyLowfb,    grPythiaNPEyHighfb);
     TGraphAsymmErrors *graMCNLO4ffb   = new TGraphAsymmErrors(grSherpa->GetN(), 
							       grX,                  graMCNLO4fYfb, 
							       grEX,                 grEX,
							       graMCNLO4fEyLowfb,    graMCNLO4fEyHighfb);
     TGraphAsymmErrors *graMCNLO4flsfb   = new TGraphAsymmErrors(grSherpa->GetN(), 
								 grX,                  graMCNLO4flsYfb, 
								 grEX,                 grEX,
								 graMCNLO4flsEyLowfb,    graMCNLO4flsEyHighfb);
     TGraphAsymmErrors *graMCNLO5fZbfb   = new TGraphAsymmErrors(grSherpa->GetN(), 
								 grX,                    graMCNLO5fZbYfb, 
								 grEX,                   grEX,
								 graMCNLO5fZbEyLowfb,    graMCNLO5fZbEyHighfb);
     TGraphAsymmErrors *graMCNLO5fZjfb   = new TGraphAsymmErrors(grSherpa->GetN(), 
								 grX,                    graMCNLO5fZjYfb, 
								 grEX,                   grEX,
								 graMCNLO5fZjEyLowfb,    graMCNLO5fZjEyHighfb);
     TGraphAsymmErrors *graMCNLO5fZb_MIfb   = new TGraphAsymmErrors(grSherpa->GetN(), 
								 grX,                    graMCNLO5fZb_MIYfb, 
								 grEX,                   grEX,
								 graMCNLO5fZb_MIEyLowfb,    graMCNLO5fZb_MIEyHighfb);
     TGraphAsymmErrors *graMCNLO5fZb_noFRfb   = new TGraphAsymmErrors(grSherpa->GetN(), 
								 grX,                    graMCNLO5fZb_noFRYfb, 
								 grEX,                   grEX,
								 graMCNLO5fZb_noFREyLowfb,    graMCNLO5fZb_noFREyHighfb);

     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);
     std::cout<<"=============================================="<<std::endl; 
     std::cout<<"====== aMC@NLO no MPI 4f UCL - v3 ============"<<std::endl;
     Double_t xsaMCNLO4f                                         = ComputeTotalCrossSection(graMCNLO4ffb);
     std::cout<<"=============================================="<<std::endl; 
     std::cout<<"====== aMC@NLO no MPI 4f Lecce - low stat ===="<<std::endl;
     Double_t xsaMCNLO4fls                                       = ComputeTotalCrossSection(graMCNLO4flsfb);
     std::cout<<"=============================================="<<std::endl; 
     std::cout<<"====== aMC@NLO no MPI 5f Zb Lecce ============"<<std::endl;
     Double_t xsaMCNLO5fZb                                       = ComputeTotalCrossSection(graMCNLO5fZbfb);
     std::cout<<"=============================================="<<std::endl; 
     std::cout<<"====== aMC@NLO no MPI 5f Zj Lecce ============"<<std::endl;
     Double_t xsaMCNLO5fZj                                       = ComputeTotalCrossSection(graMCNLO5fZjfb);
     std::cout<<"=============================================="<<std::endl; 
     std::cout<<"====== aMC@NLO no MPI 5f Zb MI Lecce ========="<<std::endl;
     Double_t xsaMCNLO5fZb_MI                                    = ComputeTotalCrossSection(graMCNLO5fZb_MIfb);
     std::cout<<"=============================================="<<std::endl; 
     std::cout<<"====== aMC@NLO no MPI 5f Zb noFR Lecce ======="<<std::endl;
     Double_t xsaMCNLO5fZb_noFR                                  = ComputeTotalCrossSection(graMCNLO5fZb_noFRfb);

     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 UCL ===================="<<std::endl;
       xsaMCNLO                                                     = ComputeTotalCrossSection(graMCNLOfb);
       std::cout<<"====== aMC@NLO NP === UCL ===================="<<std::endl;
       xsaMCNLOMPI                                                   = ComputeTotalCrossSection(graMCNLOMPIfb);
     }
     std::cout<<"=============================================="<<std::endl;
     std::cout<<"====== ALPGEN NP ============================="<<std::endl;
     Double_t xsALPGENNP=0.;
     if (grAlpgenNP!=NULL) xsALPGENNP = ComputeTotalCrossSection(grAlpgenNPfb);
     else std::cout<<"====== "<<std::endl;
     std::cout<<"====== ALPGEN NP ====no kfactor==============="<<std::endl;
     xsALPGENNP=0.;
     if (grAlpgenNP!=NULL) xsALPGENNP = ComputeTotalCrossSection(grAlpgenNPfb_nokf);
     else std::cout<<"====== "<<std::endl;
     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");
	 if (setLogY) cPlot2->SetLogy();

	 hempty->Draw();
	 
	 graMCNLOMPIfb->SetLineColor(2);
	 graMCNLOMPIfb->SetMarkerColor(2);
	 graMCNLOMPIfb->Draw("PC");
	 //	 if (showMCNLOnoMPI) 
	 //	 {
	 graMCNLOfb->SetLineColor(2);
	 graMCNLOfb->SetLineStyle(2);
	 graMCNLOfb->SetMarkerColor(2);
	 graMCNLOfb->SetMarkerStyle(24);
	 graMCNLOfb->Draw("PCSAME");

	 // UCL no MPI reprocessed with v3  
	 graMCNLO4ffb->SetLineColor(4);
	 graMCNLO4ffb->SetLineStyle(2);
	 graMCNLO4ffb->SetMarkerColor(4);
	 graMCNLO4ffb->SetMarkerStyle(20);
	 graMCNLO4ffb->Draw("PCSAME");

	 graMCNLO4flsfb->SetLineColor(2);
	 graMCNLO4flsfb->SetLineStyle(3);
	 graMCNLO4flsfb->SetMarkerColor(2);
	 graMCNLO4flsfb->SetMarkerStyle(26);
	 graMCNLO4flsfb->Draw("PCSAME");

	 graMCNLO5fZbfb->SetLineColor(7);
	 graMCNLO5fZbfb->SetLineStyle(2);
	 graMCNLO5fZbfb->SetMarkerColor(7);
	 graMCNLO5fZbfb->SetMarkerStyle(25);
	 graMCNLO5fZbfb->Draw("PCSAME");

	 graMCNLO5fZb_MIfb->SetLineColor(7);
	 graMCNLO5fZb_MIfb->SetLineStyle(1);
	 graMCNLO5fZb_MIfb->SetMarkerColor(7);
	 graMCNLO5fZb_MIfb->SetMarkerStyle(23);
	 graMCNLO5fZb_MIfb->Draw("PCSAME");


	 graMCNLO5fZb_noFRfb->SetLineColor(6);
	 graMCNLO5fZb_noFRfb->SetLineStyle(1);
	 graMCNLO5fZb_noFRfb->SetMarkerColor(6);
	 graMCNLO5fZb_noFRfb->SetMarkerStyle(25);
	 graMCNLO5fZb_noFRfb->Draw("PCSAME");

	 graMCNLO5fZjfb->SetLineColor(3);
	 graMCNLO5fZjfb->SetLineStyle(1);
	 graMCNLO5fZjfb->SetMarkerColor(3);
	 graMCNLO5fZjfb->SetMarkerStyle(21);
	 graMCNLO5fZjfb->Draw("PCSAME");
	     //	 }
	 
	 legend->Clear();
	 legend->AddEntry(graMCNLOMPIfb, "aMC@NLO authors + MPI, v<3","PEL");	
	 //if (showMCNLOnoMPI) 
	 legend->AddEntry(graMCNLOfb,     "aMC@NLO authors, v<3","PEL");	
	 legend->AddEntry(graMCNLO4ffb,   "aMC@NLO authors, v3","PEL");	
	 legend->AddEntry(graMCNLO4flsfb,   "aMC@NLO 4f LE, v3","PEL");	
	 legend->AddEntry(graMCNLO5fZbfb,   "aMC@NLO 5f Zb LE, v3","PEL");	
	 legend->AddEntry(graMCNLO5fZjfb,   "aMC@NLO 5f Zj LE, v3","PEL");	
	 legend->AddEntry(graMCNLO5fZb_MIfb,"aMC@NLO 5f Zb with MPI LE, v3","PEL");	
	 legend->AddEntry(graMCNLO5fZb_noFRfb,"aMC@NLO 5f Zb no FR. LE, v3","PEL");	
	 legend->Draw();
	 if (saveAsEps)
	   {
	     fname = "fig"+hname+"_MCNLO.eps";
	     if (setLogY) fname = "log"+fname;
	     cPlot2->SaveAs(fname.c_str());
	   }

	 TCanvas *cMCNLO_1 = new TCanvas("aMC@LNO olv vs new","aMC@LNO old vs new");
	 if (setLogY) cMCNLO_1->SetLogy();

	 hempty->Draw();

	 //	 graMCNLOMPIfb->SetLineColor(2);
	 //	 graMCNLOMPIfb->SetMarkerColor(2);
	 graMCNLOMPIfb->Draw("PC");
	 //	 if (showMCNLOnoMPI) 
	 //	 {
	 //	 graMCNLOfb->SetLineColor(2);
	 //	 graMCNLOfb->SetLineStyle(2);
	 //	 graMCNLOfb->SetMarkerColor(2);
	 //	 graMCNLOfb->SetMarkerStyle(24);
	 graMCNLOfb->Draw("PCSAME");

	 // UCL no MPI reprocessed with v3  
	 //	 graMCNLO4ffb->SetLineColor(4);
	 //	 graMCNLO4ffb->SetLineStyle(2);
	 //	 graMCNLO4ffb->SetMarkerColor(4);
	 //	 graMCNLO4ffb->SetMarkerStyle(20);
	 graMCNLO4ffb->Draw("PCSAME");

	 //	 graMCNLO4flsfb->SetLineColor(2);
	 //	 graMCNLO4flsfb->SetLineStyle(3);
	 //	 graMCNLO4flsfb->SetMarkerColor(2);
	 //	 graMCNLO4flsfb->SetMarkerStyle(26);
	 graMCNLO4flsfb->Draw("PCSAME");

	 legend->Draw();
	 if (saveAsEps)
	   {
	     fname = "fig"+hname+"_MCNLO_all4f.eps";
	     if (setLogY) fname = "log"+fname;
	     cMCNLO_1->SaveAs(fname.c_str());
	   }


	 TCanvas *cMCNLO_2 = new TCanvas("aMC@LNO 4f vs 5f no MPI","aMC@LNO 4f vs 5f no MPI");
	 if (setLogY) cMCNLO_2->SetLogy();
	 hempty->Draw();

	 graMCNLOfb->Draw("PCSAME");

	 // UCL no MPI reprocessed with v3  
	 graMCNLO4ffb->Draw("PCSAME");

	 graMCNLO4flsfb->Draw("PCSAME");

	 graMCNLO5fZbfb->Draw("PCSAME");

	 legend->Draw();
	 if (saveAsEps)
	   {
	     fname = "fig"+hname+"_MCNLO_4fvs5f_noMPI.eps";
	     if (setLogY) fname = "log"+fname;
	     cMCNLO_2->SaveAs(fname.c_str());
	   }


	 TCanvas *cMCNLO_3 = new TCanvas("aMC@LNO 4f vs 5f w MPI","aMC@LNO 4f vs 5f w MPI");
	 if (setLogY) cMCNLO_3->SetLogy();
	 hempty->Draw();

	 graMCNLOMPIfb->Draw("PC");
	 graMCNLO5fZb_MIfb->Draw("PCSAME");

	 legend->Draw();
	 if (saveAsEps)
	   {
	     fname = "fig"+hname+"_MCNLO_4fvs5f_wMPI.eps";
	     if (setLogY) fname = "log"+fname;
	     cMCNLO_3->SaveAs(fname.c_str());
	   }



	 TCanvas *cMCNLO_4 = new TCanvas("aMC@LNO 5f Zb vs 5f Zj wo MPI","aMC@LNO 5f Zb vs 5f Zj wo MPI");
	 if (setLogY) cMCNLO_4->SetLogy();
	 hempty->Draw();

	 graMCNLOfb->Draw("PC");
	 graMCNLO4ffb->Draw("PCSAME");
	 graMCNLO5fZbfb->Draw("PCSAME");
	 graMCNLO5fZjfb->Draw("PCSAME");

	 legend->Draw();
	 if (saveAsEps)
	   {
	     fname = "fig"+hname+"_MCNLO_5f_ZbvsZj_wMPI.eps";
	     if (setLogY) fname = "log"+fname;
	     cMCNLO_3->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();

     grMCFMNPsta->SetFillColor(colfMCFMstat);
     grMCFM8NPsta->SetFillColor(colfMCFMstat);
     grMCFM23NPsta->SetFillColor(colfMCFMstat);
     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 (grAlpgenNP!=NULL)
     {
	 grAlpgenNPfb->SetLineStyle(1);
	 grAlpgenNPfb->SetLineColor(colAlpgen);
	 grAlpgenNPfb->SetMarkerColor(colAlpgen);
	 grAlpgenNPfb->Draw("PCSAME");
	 grAlpgenNPfb_nokf->SetLineStyle(2);
	 grAlpgenNPfb_nokf->SetLineColor(colAlpgen);
	 grAlpgenNPfb_nokf->SetMarkerColor(colAlpgen);
	 grAlpgenNPfb_nokf->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 (grAlpgenNP!=NULL) {
	 legend->AddEntry(grAlpgenNPfb,"Alpgen","PEL");
	 legend->AddEntry(grAlpgenNPfb_nokf,"Alpgen no k-factor","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;



     TCanvas *cPlotSA = new TCanvas("SherpaAlpgenAmcnlo","SherpaAlpgenAmcnlo", 500,700);
     //////
     //////
     // 
     Double_t pBottSize = 0.45;
     Double_t pTopSize  = 1.-pBottSize;
     Double_t BotMarg = 0.2;
     Double_t TopMarg = 0.05;
     Double_t pBottSizeEff = pBottSize-BotMarg;
     Double_t pTopSizeEff  = pTopSize-TopMarg;
     TPad *ct_1 = new TPad("ct_1", "ctnewpad", 0.0,pBottSize,1.0,1.0);
     ct_1->SetTopMargin(   TopMarg );
     ct_1->SetBottomMargin(0.0);
     ct_1->SetFillStyle(c1->GetFillStyle());
     ct_1->Draw();
     TPad *ct_2 = new TPad("ct_2", "ctnewpad2",0.0,0.0,1.0,pBottSize);
     ct_2->SetFillStyle(c1->GetFillStyle());
     ct_2->SetTopMargin   (0.0);
     ct_2->SetBottomMargin(BotMarg);
     ct_2->Draw();

     if (setLogY) {
	 ct_1->SetLogy();
	 //ct_2->SetLogx();
     }     

     ct_1->cd();
     //hempty->Draw();
     TH2F* hemptyCopy = hempty->Clone();
     double sizeTopPlot = 0.05;
     hemptyCopy->GetXaxis()->SetLabelSize(0.);
     hemptyCopy->GetXaxis()->SetTitleSize(0.);
     hemptyCopy->GetYaxis()->SetLabelSize(sizeTopPlot);
     hemptyCopy->GetYaxis()->SetTitleSize(sizeTopPlot);
     hemptyCopy->GetYaxis()->SetTitleOffset(1.2);
     double sizeBottomPlot = sizeTopPlot*pTopSize/pBottSize;
     hemptyRat->GetYaxis()->SetLabelSize(sizeBottomPlot);
     hemptyRat->GetYaxis()->SetTitleSize(sizeBottomPlot);
     hemptyRat->GetXaxis()->SetLabelSize(sizeBottomPlot);
     hemptyRat->GetXaxis()->SetTitleSize(sizeBottomPlot);
     hemptyRat->GetYaxis()->SetTitleOffset(1.);
     hemptyCopy->Draw();
     //////
     //////
     //grSherpafb->Draw("PC");
     //grPythiafb->Draw("PCSAME");
     grSherpaNPfb->Draw("PCSAME");
     //grPythiaNPfb->Draw("PCSAME");
     if (grAlpgenNP!=NULL)
     {
	 grAlpgenNPfb_nokf->Draw("PCSAME");
     }     
     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO)) {
       graMCNLOMPIfb->Draw("PCSAME");
//        if (showMCNLOnoMPI) 
//        {
// 	   graMCNLOfb->Draw("PCSAME");
//        }
     }
       graMCNLO5fZb_MIfb->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 (grAlpgenNP!=NULL) legend->AddEntry(grAlpgenNPfb_nokf,"Alpgen no k-factor","PEL");
     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))  
     {
	 //if (showMCNLOnoMPI) legend->AddEntry(graMCNLOfb,   "aMC@NLO no MPI","PEL");	
       legend->AddEntry(graMCNLOMPIfb, "aMC@NLO + MPI","PEL");	
     }
     legend->AddEntry(graMCNLO5fZb_MIfb, "aMC@NLO 5f Zb + MPI","PEL");	
     legend->Draw();
     
     // 
     // 
     ct_2->cd();
     hemptyRat->Draw();
     TGraphAsymmErrors* mlrat ;
     TGraphAsymmErrors* nlorat ;
     nlorat = getGraphRatio(graMCNLOMPIfb,graMCNLO5fZb_MIfb);
     if (grAlpgenNP!=NULL) {
	 mlrat  = getGraphRatio(grAlpgenNPfb_nokf,grSherpaNPfb);
	 mlrat->Draw("PESAME");
     }
     nlorat->Draw("PESAME");
     legendRATD->Clear();
      if (grAlpgenNP!=NULL) legendRATD->AddEntry( mlrat,"multi-leg 5f/4f", "PEL");
     legendRATD->AddEntry(nlorat,"NLO       5f/4f", "PEL");
     legendRATD->Draw();

     if (saveAsEps)
       {
	 fname = "fig"+hname+"_aMCNLO45f_Sherpa_Alpgen.eps";
	 if (setLogY) fname = "log"+fname;
	 cPlotSA->SaveAs(fname.c_str());
       }

     TCanvas *cPlotSPA = new TCanvas("SherpaPyhiaAlpgenAmcnloMCFM","SherpaPythiaAlpgenAmcnloMCFM");
     if (setLogY) cPlotSPA->SetLogy();
     
     hempty->Draw();

     //grSherpafb->Draw("PC");
     //grPythiafb->Draw("PCSAME");
     grSherpaNPfb->Draw("PCSAME");
     grPythiaNPfb->Draw("PCSAME");
     if (grAlpgenNP!=NULL)
     {
	 grAlpgenNPfb_nokf->Draw("PCSAME");
	 grAlpgenNPfb->Draw("PCSAME");
     }     
     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO)) {
       graMCNLOMPIfb->Draw("PCSAME");
//        if (showMCNLOnoMPI) 
//        {
// 	   graMCNLOfb->Draw("PCSAME");
//        }
     }
       graMCNLO5fZb_MIfb->Draw("PCSAME");
       grMCFM8NPsta_as_pdf_sca_np->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 (grAlpgenNP!=NULL) {
	 legend->AddEntry(grAlpgenNPfb,     "Alpgen with k-factor","PEL");
	 legend->AddEntry(grAlpgenNPfb_nokf,"Alpgen no k-factor","PEL");
     }
     if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))  
     {
	 //if (showMCNLOnoMPI) legend->AddEntry(graMCNLOfb,   "aMC@NLO no MPI","PEL");	
       legend->AddEntry(graMCNLOMPIfb, "aMC@NLO authors + MPI","PEL");	
     }
     legend->AddEntry(graMCNLO5fZb_MIfb, "aMC@NLO 5f Zb + MPI","PEL");	
     legend->AddEntry(grMCFM8NPsta_as_pdf_sca_np, "MCFM mstw8 - total err", "PEL");

     legend->Draw();
     
     if (saveAsEps)
       {
	 fname = "fig"+hname+"_aMCNLO45f_Sherpa_Py_Alpgen.eps";
	 if (setLogY) fname = "log"+fname;
	 cPlotSPA->SaveAs(fname.c_str());
       }


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


     ///////+++++++++++++++++++++++++++++++++++++++ try with ratio 
     //if (cPlotNP) delete cPlotNP;
     TCanvas *cPlotNP_2canvas = new TCanvas("MCFMwNP_withRatio","MCFMwNP_withRatio", 500, 700);
     Double_t pBottSize = 0.45;
     Double_t pTopSize  = 1.-pBottSize;
     Double_t BotMarg = 0.2;
     Double_t TopMarg = 0.05;
     Double_t pBottSizeEff = pBottSize-BotMarg;
     Double_t pTopSizeEff  = pTopSize-TopMarg;
     TPad *c1_1 = new TPad("c1_1", "newpad", 0.0,pBottSize,1.0,1.0);
     c1_1->SetTopMargin(   TopMarg );
     c1_1->SetBottomMargin(0.0);
     c1_1->SetFillStyle(c1->GetFillStyle());
     c1_1->Draw();
     TPad *c1_2 = new TPad("c1_2", "newpad2",0.0,0.0,1.0,pBottSize);
     c1_2->SetFillStyle(c1->GetFillStyle());
     c1_2->SetTopMargin   (0.0);
     c1_2->SetBottomMargin(BotMarg);
     c1_2->Draw();

     if (setLogY) c1_1->SetLogy();
     c1_1->cd();
     //hempty->Draw();
     //TH2F* hemptyCopy = hempty->Clone();
     double sizeTopPlot = 0.05;
     hemptyCopy->GetXaxis()->SetLabelSize(0.);
     hemptyCopy->GetXaxis()->SetTitleSize(0.);
     hemptyCopy->GetYaxis()->SetLabelSize(sizeTopPlot);
     hemptyCopy->GetYaxis()->SetTitleSize(sizeTopPlot);
     hemptyCopy->GetYaxis()->SetTitleOffset(1.2);
     double sizeBottomPlot = sizeTopPlot*pTopSize/pBottSize;
     hemptyRat->GetYaxis()->SetLabelSize(sizeBottomPlot);
     hemptyRat->GetYaxis()->SetTitleSize(sizeBottomPlot);
     hemptyRat->GetXaxis()->SetLabelSize(sizeBottomPlot);
     hemptyRat->GetXaxis()->SetTitleSize(sizeBottomPlot);
     hemptyRat->GetYaxis()->SetTitleOffset(1.);
     hemptyCopy->Draw();



     //return;
     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("PE2SAME");
     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();
     
     //     return;
     // Bottom 
     //if (setLogY && ("j_bjet_ZdR620_binned"==hname)) c1_2->SetLogy();
     c1_2->cd();

     hemptyRat->Draw();
     TGraphAsymmErrors* rat0=getGraphRatio(grMCFMNPsta_as_pdf_sca_np,grMCFMNPsta_as_pdf_sca_np);
     TGraphAsymmErrors* rat1=getGraphRatio(grMCFMNPsta_as_pdf_sca_np,grMCFMNPsta);
     TGraphAsymmErrors* rat2=getGraphRatio(grMCFMNPsta_as_pdf_sca_np,graMCNLOMPIfb);
     TGraphAsymmErrors* rat3=getGraphRatio(grMCFMNPsta_as_pdf_sca_np,grPythiaNPfb);
     TGraphAsymmErrors* rat4=getGraphRatio(grMCFMNPsta_as_pdf_sca_np,grSherpaNPfb);

     rat0->Draw("E2");
     rat1->Draw("E2SAME");
     rat2->Draw("EPSAME");
     rat3->Draw("EPSAME");
     rat4->Draw("EPSAME");


     if (saveAsEps)
       {
	 fname = "fig"+hname+"_MCFM_aMCNLO_Sherpa_Pythia_WITHRATIO.eps";
	 if (setLogY) fname = "log"+fname;
	 cPlotNP_2canvas->SaveAs(fname.c_str());
       }
     //     return;

     
     //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("PE2SAME");
     grMCFM8NPsta->Draw("PESAME");

     // CT10 central value with NP 
     grMCFMNPsta->Draw("PESAME");
     grMCFM23NPsta->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(grMCFMNPsta            ,   "MCFM\_nnpdf23 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("PE2SAME");
     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);
	 grdataMbbb->SetMarkerColor(1);
	 grdataNbbb->SetMarkerStyle(24);
	 grdataNbbb->SetMarkerSize(1.5);
	 grdataMbbb->SetMarkerStyle(26);
	 grdataMbbb->SetMarkerSize(1.5);


// 	 grdataNBayesSyst->SetLineColor(2);
// 	 grdataNBayesSyst->SetLineWidth(2);
// 	 grdataMBayesSyst->SetLineColor(2);
// 	 grdataMBayesSyst->SetLineWidth(2);

	 grdataMBayesSyst->SetMarkerColor(1);
	 grdataNBayesSyst->SetMarkerStyle(20);
	 grdataMBayesSyst->SetMarkerStyle(23);
	 grdataNBayesSyst->SetMarkerSize(0.01);
	 grdataMBayesSyst->SetMarkerSize(0.01);
	 grdataNBayesTot->SetMarkerColor(1);
	 grdataMBayesTot->SetMarkerColor(1);
	 grdataNBayesTot->SetMarkerStyle(20);
	 grdataMBayesTot->SetMarkerStyle(23);
	 grdataNBayesTot->SetMarkerSize(0.01);
	 grdataMBayesTot->SetMarkerSize(0.01);
    

	 grdataMbbb->Draw("PESAME");
	 grdataMBayes->Draw("PESAME");     
	 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");
	     legend1->AddEntry(grdataNbbb,               "Nicola bxb Alpgen-unf.","PEL");
	     legend1->AddEntry(grdataNBayes,             "Nicola Bayes Alpgen-unf.","PEL");
	   }
	 legend1->Draw();
     


	 if (saveAsEps)
	   {
	       if (MBayes) fname = "fig"+hname+"_MCFM_w_wo_NP_data.eps";
	       else if (NBayes) fname = "fig"+hname+"_MCFM_w_wo_NP_dataNic.eps";
	     if (setLogY) fname = "log"+fname;
	     cPlotNP_01data->SaveAs(fname.c_str());
	   }



	 TCanvas *cDandR = new TCanvas("DataAndRatio","DataAndRatio", 500, 700);
	 TPad *c1DandR = new TPad("c1_1data", "newpad", 0.0,pBottSize,1.0,1.0);
	 c1DandR->SetTopMargin(   TopMarg );
	 c1DandR->SetBottomMargin(0.0);
	 c1DandR->SetFillStyle(c1->GetFillStyle());
	 c1DandR->Draw();
	 TPad *c2DandR = new TPad("c1_2data", "newpad2",0.0,0.0,1.0,pBottSize);
	 c2DandR->SetFillStyle(c1->GetFillStyle());
	 c2DandR->SetTopMargin   (0.0);
	 c2DandR->SetBottomMargin(BotMarg);
	 c2DandR->Draw();

	 if (setLogY) c1DandR->SetLogy();
	 c1DandR->cd();
	 hemptyCopy->Draw();


// 	 grMCFMNPsta_as_pdf_sca_np->Draw("E2");
// 	 grMCFM8NPsta_as_pdf_sca_np->Draw("E2SAME");
// 	 grMCFM23NPsta_as_pdf_sca_np->Draw("E2SAME");
	 grMCFM8NPsta_as_pdf_sca_np->SetLineColor(4);	 
	 grMCFMNPsta_as_pdf_sca_np->SetLineColor(3);	 
	 grMCFM23NPsta_as_pdf_sca_np->SetLineColor(7);	 

	 grMCFM8NPsta_as_pdf_sca_np ->SetLineWidth(2);	 
	 grMCFMNPsta_as_pdf_sca_np  ->SetLineWidth(2);	 
	 grMCFM23NPsta_as_pdf_sca_np->SetLineWidth(2);	 
	 //
	 grMCFM8NPsta_as_pdf_sca_np->SetMarkerColor(4);	 
	 grMCFMNPsta_as_pdf_sca_np->SetMarkerColor(3);	 
	 grMCFM23NPsta_as_pdf_sca_np->SetMarkerColor(7);	 
	 //
	 grMCFM8NPsta_as_pdf_sca_np->SetMarkerSize(0.01);	 
	 grMCFMNPsta_as_pdf_sca_np->SetMarkerSize(0.01);	 
	 grMCFM23NPsta_as_pdf_sca_np->SetMarkerSize(0.01);	 
	 //
	 //grMCFMNPsta_as_pdf_sca_np  ->Draw("CX");
	 grMCFM8NPsta_as_pdf_sca_np ->Draw("CX");
	 grSherpaNPfb ->SetLineWidth(3);
	 grSherpaNPfb ->SetMarkerSize(0.01);
	 grSherpaNPfb ->Draw("CX");
	 //
	 if (grAlpgenNP!=NULL)
	 {
	     grAlpgenNPfb ->SetLineWidth(3);
	     grAlpgenNPfb ->SetMarkerSize(0.01);
	     grAlpgenNPfb ->Draw("CX");
	 }
	 
	 //grMCFM8NPsta_as_pdf_sca_np ->Draw("E4");
	 //grMCFM23NPsta_as_pdf_sca_np->Draw("CX");

	 if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))   
	   {     
	     /// test ---------------{
	     graMCNLOMPIfb->SetMarkerSize(0.01);
	     graMCNLOMPIfb->SetLineWidth(2);
	     graMCNLOMPIfb->SetLineStyle(1);
	     graMCNLOMPIfb->SetLineColor(2);
	     graMCNLOMPIfb->SetMarkerColor(2);
	     graMCNLOMPIfb->TGraphPainter::Draw("CX");
	     /// end test -----------}
	     /// test ---------------{
	     if (showMCNLOnoMPI)
	       {
		 graMCNLOfb->SetLineColor(2);
		 graMCNLOfb->SetLineStyle(2);
		 graMCNLOfb->SetMarkerColor(2);
		 graMCNLOfb->SetMarkerStyle(24);
		 graMCNLOfb->Draw("PCSAME");
	       }
	     /// end test -----------}
	   }

	 grdataNbbb  ->SetFillColor(5);
	 grdataNBayes->SetFillColor(5);
	 grdataMbbb  ->SetFillColor(5);
	 grdataMBayes->SetFillColor(5);
	 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);
    

	 if (Mbbb)   grdataMbbb->Draw("PESAME");
	 if (MBayes) {
	     grdataMBayesTot->Draw("PESAME");     
	     //grdataMBayesSyst->Draw("PESAME");     
	 }
	 if (Nbbb)   grdataNbbb->Draw("PESAME");
	 if (NBayes) {
	     grdataNBayesTot->Draw("PESAME");
	     //grdataNBayesSyst->Draw("PESAME");
	 }
	 
	 legendDMC->Clear();
 	 legendDMC->AddEntry(grMCFM8NPsta_as_pdf_sca_np,"MCFM (mstw2008) x NP","L");
	 //legendDMC->AddEntry(grMCFMNPsta_as_pdf_sca_np,"MCFM (ct10) x NP","L");
	 //legendDMC->AddEntry(grMCFM23NPsta_as_pdf_sca_np,"MCFM (nnpdf23) x NP","L");
	 if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))  
	   {
	     legendDMC->AddEntry(graMCNLOMPIfb,            "aMC@NLO","L");
	     if (showMCNLOnoMPI) legendDMC->AddEntry(graMCNLOfb,               "aMC@NLO no MPI","L");
	   }
	 legendDMC->AddEntry(grSherpaNPfb,"SHERPA","L");
	 if (grAlpgenNP!=NULL) legendDMC->AddEntry(grAlpgenNPfb,"Alpgen","L");
	 
	 if (MBayes) legendDMC->AddEntry(grdataMBayes,             "ATLAS data, 4.6 fb^{-1}","PEL");
	 if (NBayes) legendDMC->AddEntry(grdataNBayes,             "ATLAS data, 4.6 fb^{-1}","PEL");
	 if (Mbbb)   legendDMC->AddEntry(grdataMbbb,               "ATLAS data, 4.6 fb^{-1}","PEL");
	 if (Nbbb)   legendDMC->AddEntry(grdataNbbb,               "ATLAS data, 4.6 fb^{-1}","PEL");	 
	 legendDMC->Draw();
	 if (hname=="j_bjet_ZdPhiF20_binned") 
	 {
	     ATLAS_LABEL(0.2,0.55);
	     myText(0.33,0.55,1,"Internal");
	     myText(0.65,0.15,1,"Z+b selection");	     
	 }
	 else if (hname=="j_bjet_ZdR620_binned") 
	 {
	     ATLAS_LABEL(0.2,0.55);
	     myText(0.33,0.55,1,"Internal");
	     myText(0.65,0.25,1,"Z+b selection");	     
	 }
	 else 
	 {
	     ATLAS_LABEL(0.2,0.87);
	     myText(0.33,0.87,1,"Internal");
	     myText(0.6,0.5,1,"Z+b selection");
	 }
	 

	 c2DandR->cd();
	 /*
	 Double_t ey20up[100];
	 Double_t ey20dw[100];
	 for (int j=0; j<grdataMBayes->GetN(); ++j) {
	     ey20up[j]=0.15*grdataMBayes->GetY()[j];
	     ey20dw[j]=0.15*grdataMBayes->GetY()[j];
	 }
	 
	 Double_t eylN[100];
	 Double_t eyhN[100];
	 Double_t eylM[100];
	 Double_t eyhM[100];
	 for (int j=0; j<grdataMBayes->GetN(); ++j) 
	 {
	     eylM[j]=sqrt( edataMBayesSdw[j]*edataMBayesSdw[j] + (grdataMBayes->GetEYlow()[j] *grdataMBayes->GetEYlow()[j])  );
	     eyhM[j]=sqrt( edataMBayesSup[j]*edataMBayesSup[j] + (grdataMBayes->GetEYhigh()[j]*grdataMBayes->GetEYhigh()[j]) );

	     eylN[j]=sqrt( edataNBayesSdw[j]*edataNBayesSdw[j] + (grdataNBayes->GetEYlow()[j] *grdataNBayes->GetEYlow()[j])  );
	     eyhN[j]=sqrt( edataNBayesSup[j]*edataNBayesSup[j] + (grdataNBayes->GetEYhigh()[j]*grdataNBayes->GetEYhigh()[j]) );

	 }
	 */
	 
	 // syst only - assuming 15% 
	 TGraphAsymmErrors* grdataSyst;
	 // stat + syst		      
	 TGraphAsymmErrors* grdataStat;

	 if (MBayes) {
// 	   grdataSyst = new TGraphAsymmErrors(grdataMBayes->GetN(),
// 					      grdataMBayes->GetX(),     grdataMBayes->GetY(),
// 					      grdataMBayes->GetEXlow(), grdataMBayes->GetEXhigh(),
// 					      edataMBayesSdw,           edataMBayesSup);

// 	   grdataStat = new TGraphAsymmErrors(grdataMBayes->GetN(),
// 					      grdataMBayes->GetX(),     grdataMBayes->GetY(),
// 					      grdataMBayes->GetEXlow(), grdataMBayes->GetEXhigh(),
// 					      eylM, eyhM);
	     grdataSyst = grdataMBayesSyst;	     
	     grdataStat = grdataMBayesTot;	     
	 }
	 if (NBayes) {
// 	   grdataSyst = new TGraphAsymmErrors(grdataNBayes->GetN(),
// 					      grdataNBayes->GetX(),     grdataNBayes->GetY(),
// 					      grdataNBayes->GetEXlow(), grdataNBayes->GetEXhigh(),
// 					      edataNBayesSdw,           edataNBayesSup);

// 	   grdataStat = new TGraphAsymmErrors(grdataNBayes->GetN(),
// 					      grdataNBayes->GetX(),     grdataNBayes->GetY(),
// 					      grdataNBayes->GetEXlow(), grdataNBayes->GetEXhigh(),
// 					      eylN, eyhN);
	     grdataSyst = grdataNBayesSyst;	     
	     grdataStat = grdataNBayesTot;	     
	 }
	 
	 grdataSyst->SetFillColor(1);
	 grdataSyst->SetFillStyle(3005);
	 grdataStat->SetFillColor(5);
	 

     //     if (setLogY) c1_1->SetLogy();
	 hemptyRat->Draw();
// 	 grMCFMNPsta_as_pdf_sca_np->SetMarkerColor(colmMCFM);
// 	 grMCFMNPsta_as_pdf_sca_np->SetLineColor(colmMCFM);
// 	 grMCFM8NPsta_as_pdf_sca_np->SetMarkerColor(colmMCFM8);
// 	 grMCFM8NPsta_as_pdf_sca_np->SetLineColor(colmMCFM8);
// 	 grMCFM23NPsta_as_pdf_sca_np->SetMarkerColor(colmMCFM23);
// 	 grMCFM23NPsta_as_pdf_sca_np->SetLineColor(colmMCFM23);
	 


	 TGraphAsymmErrors* rat0D ;
	 TGraphAsymmErrors* rat1D ;
	 TGraphAsymmErrors* rat2D ;
	 TGraphAsymmErrors* rat3D ;
	 TGraphAsymmErrors* rat4D ;
	 TGraphAsymmErrors* rat5D ;
	 TGraphAsymmErrors* rat6D ;
	 TGraphAsymmErrors* rat7D ;

	 if (MBayes) {
	   rat0D = getGraphRatio(grdataMBayes,grdataStat);
	   rat1D = getGraphRatio(grdataMBayes,grdataSyst);
	   rat2D = getGraphRatio(grdataMBayes,grMCFMNPsta_as_pdf_sca_np);
	   rat3D = getGraphRatio(grdataMBayes,grMCFM8NPsta_as_pdf_sca_np);
	   rat4D = getGraphRatio(grdataMBayes,grMCFM23NPsta_as_pdf_sca_np);
	   rat5D = getGraphRatio(grdataMBayes,graMCNLOMPIfb);
	   rat6D = getGraphRatio(grdataMBayes,grSherpaNPfb);
	   if (grAlpgenNP!=NULL) rat7D = getGraphRatio(grdataMBayes,grAlpgenNPfb);
	 }
	 if (NBayes) {
	   rat0D = getGraphRatio(grdataNBayes,grdataStat);
	   rat1D = getGraphRatio(grdataNBayes,grdataSyst);
	   rat2D = getGraphRatio(grdataNBayes,grMCFMNPsta_as_pdf_sca_np);
	   rat3D = getGraphRatio(grdataNBayes,grMCFM8NPsta_as_pdf_sca_np);
	   rat4D = getGraphRatio(grdataNBayes,grMCFM23NPsta_as_pdf_sca_np);
	   rat5D = getGraphRatio(grdataNBayes,graMCNLOMPIfb);
	   rat6D = getGraphRatio(grdataNBayes,grSherpaNPfb);
	   if (grAlpgenNP!=NULL) rat7D = getGraphRatio(grdataNBayes,grAlpgenNPfb);
	 }

	 rat0D->Draw("E2");
	 rat1D->Draw("E2SAME");
	 //rat2D->Draw("EPSAME");
	 rat5D->Draw("EPSAME");
	 rat6D->Draw("EPSAME");
	 rat7D->Draw("EPSAME");
	 rat3D->Draw("EPSAME");
	 //rat4D->Draw("EPSAME");
	 legendRATD->Clear();
 	 legendRATD->AddEntry(rat0D,"Data stat+syst err.","F");
 	 legendRATD->AddEntry(rat1D,"Data syst err.","F");
	 legendRATS->Clear();
	 legendRATS->AddEntry(rat3D,"MCFM(mstw8)xNP, stat+syst err.","L");
	 //legendRATS->AddEntry(rat2D,"MCFM(ct10)xNP, stat+syst err. ","L");
	 //legendRATS->AddEntry(rat4D,"MCFM(nnpdf23)xNP, stat+syst err. ","L");
	 legendRATS->AddEntry(rat5D,"aMC@NLO","L");
	 legendRATS->AddEntry(rat6D,"SHERPA","L");
	 if (grAlpgenNP!=NULL) legendRATS->AddEntry(rat7D,"Alpgen","L");
	 legendRATD->Draw();
	 legendRATS->Draw();
	 
	 if (saveAsEps)
	   {
	       if (MBayes) fname = "fig"+hname+"_MCFM_w_wo_NP_data_WITHRATIO.eps";
	       else if (NBayes) fname = "fig"+hname+"_MCFM_w_wo_NP_dataNic_WITHRATIO.eps"; 
	     if (setLogY) fname = "log"+fname;
	     cDandR->SaveAs(fname.c_str());
	   }

	 //	 return;

	 ///////////////////////////// reminder {
     // Double_t pBottSize = 0.45;
     // Double_t pTopSize  = 1.-pBottSize;
     // Double_t BotMarg = 0.2;
     // Double_t TopMarg = 0.05;
     // Double_t pBottSizeEff = pBottSize-BotMarg;
     // Double_t pTopSizeEff  = pTopSize-TopMarg;
     // double sizeTopPlot = 0.05;
     // hemptyCopy->GetXaxis()->SetLabelSize(0.);
     // hemptyCopy->GetXaxis()->SetTitleSize(0.);
     // hemptyCopy->GetYaxis()->SetLabelSize(sizeTopPlot);
     // hemptyCopy->GetYaxis()->SetTitleSize(sizeTopPlot);
     // hemptyCopy->GetYaxis()->SetTitleOffset(1.2);
     // double sizeBottomPlot = sizeTopPlot*pTopSize/pBottSize;
     // hemptyRat->GetYaxis()->SetLabelSize(sizeBottomPlot);
     // hemptyRat->GetYaxis()->SetTitleSize(sizeBottomPlot);
     // hemptyRat->GetXaxis()->SetLabelSize(sizeBottomPlot);
     // hemptyRat->GetXaxis()->SetTitleSize(sizeBottomPlot);
     // hemptyRat->GetYaxis()->SetTitleOffset(1.);
	 ///////////////////////////// reminder }


	 TCanvas *cDandRPDF = new TCanvas("DataAndRatioPDF","DataAndRatioPDF", 500, 700);
	 //
	 TPad *c1DandRPDF = new TPad("c1_1dataPDF", "newpadPDF", 0.0,0.6,1.0,1.0);
	 c1DandRPDF->SetTopMargin(   TopMarg );
	 c1DandRPDF->SetBottomMargin(0.0);
	 c1DandRPDF->SetFillStyle(c1->GetFillStyle());
	 c1DandRPDF->Draw();
	 TPad *c2DandRPDF = new TPad("c1_2dataPDF", "newpad2PDF",0.0,0.31,1.0,0.6);
	 c2DandRPDF->SetFillStyle(c1->GetFillStyle());
	 c2DandRPDF->SetTopMargin   (0.0);
	 c2DandRPDF->SetBottomMargin(0.0);
	 c2DandRPDF->Draw();
	 TPad *c3DandRPDF = new TPad("c1_3dataPDF", "newpad3PDF",0.0,0.0,1.0,0.31);
	 c3DandRPDF->SetFillStyle(c1->GetFillStyle());
	 c3DandRPDF->SetTopMargin   (0.0);
	 c3DandRPDF->SetBottomMargin(BotMarg);
	 c3DandRPDF->Draw();

	 if (setLogY) c1DandRPDF->SetLogy();
	 c1DandRPDF->cd();
	 hemptyCopy->Draw();


// 	 grMCFMNPsta_as_pdf_sca_np->Draw("E2");
// 	 grMCFM8NPsta_as_pdf_sca_np->Draw("E2SAME");
// 	 grMCFM23NPsta_as_pdf_sca_np->Draw("E2SAME");
	 grMCFM8NPsta_as_pdf_sca_np->SetLineColor(4);	 
	 grMCFMNPsta_as_pdf_sca_np->SetLineColor(3);	 
	 grMCFM23NPsta_as_pdf_sca_np->SetLineColor(7);	 

	 grMCFM8NPsta_as_pdf_sca_np ->SetLineWidth(2);	 
	 grMCFMNPsta_as_pdf_sca_np  ->SetLineWidth(2);	 
	 grMCFM23NPsta_as_pdf_sca_np->SetLineWidth(2);	 
	 //
	 grMCFM8NPsta_as_pdf_sca_np->SetMarkerColor(4);	 
	 grMCFMNPsta_as_pdf_sca_np->SetMarkerColor(3);	 
	 grMCFM23NPsta_as_pdf_sca_np->SetMarkerColor(7);	 
	 //
	 grMCFM8NPsta_as_pdf_sca_np->SetMarkerSize(0.01);	 
	 grMCFMNPsta_as_pdf_sca_np->SetMarkerSize(0.01);	 
	 grMCFM23NPsta_as_pdf_sca_np->SetMarkerSize(0.01);	 
	 //
	 grMCFMNPsta_as_pdf_sca_np  ->Draw("CX");
	 grMCFM8NPsta_as_pdf_sca_np ->Draw("CX");
	 grSherpaNPfb ->SetLineWidth(3);
	 grSherpaNPfb ->SetMarkerSize(0.01);
	 //grSherpaNPfb ->Draw("CX");
	 grMCFM23NPsta_as_pdf_sca_np->Draw("CX");

	 if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))   
	   {     
	     /// test ---------------{
	     graMCNLOMPIfb->SetMarkerSize(0.01);
	     graMCNLOMPIfb->SetLineWidth(2);
	     graMCNLOMPIfb->SetLineStyle(1);
	     graMCNLOMPIfb->SetLineColor(2);
	     graMCNLOMPIfb->SetMarkerColor(2);
	     //graMCNLOMPIfb->TGraphPainter::Draw("CX");
	     /// end test -----------}
	     /// test ---------------{
	     if (showMCNLOnoMPI)
	       {
		 graMCNLOfb->SetLineColor(2);
		 graMCNLOfb->SetLineStyle(2);
		 graMCNLOfb->SetMarkerColor(2);
		 graMCNLOfb->SetMarkerStyle(24);
		 //graMCNLOfb->Draw("PCSAME");
	       }
	     /// end test -----------}
	   }

	 grdataNbbb  ->SetFillColor(5);
	 grdataNBayes->SetFillColor(5);
	 grdataMbbb  ->SetFillColor(5);
	 grdataMBayes->SetFillColor(5);
	 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);
    

	 if (Mbbb)   grdataMbbb->Draw("PESAME");
	 if (MBayes) {
	     grdataMBayesTot->Draw("PESAME");     
	     //grdataMBayesSyst->Draw("PESAME");     
	 }
	 if (Nbbb)   grdataNbbb->Draw("PESAME");
	 if (NBayes) {
	     grdataNBayesTot->Draw("PESAME");
	     //grdataNBayesSyst->Draw("PESAME");
	 }
	 
	 legendDMC->Clear();
 	 legendDMC->AddEntry(grMCFM8NPsta_as_pdf_sca_np,"MCFM (mstw2008) x NP","L");
	 legendDMC->AddEntry(grMCFMNPsta_as_pdf_sca_np,"MCFM (ct10) x NP","L");
	 legendDMC->AddEntry(grMCFM23NPsta_as_pdf_sca_np,"MCFM (nnpdf23) x NP","L");
	 //if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))  
	 //  {
	     //legendDMC->AddEntry(graMCNLOMPIfb,            "aMC@NLO","L");
	     //if (showMCNLOnoMPI) legendDMC->AddEntry(graMCNLOfb,               "aMC@NLO no MPI","L");
	 //  }
	 //legendDMC->AddEntry(grSherpaNPfb,"SHERPA","L");
	 if (MBayes) legendDMC->AddEntry(grdataMBayes,             "ATLAS data, 4.6 fb^{-1}","PEL");
	 if (NBayes) legendDMC->AddEntry(grdataNBayes,             "ATLAS data, 4.6 fb^{-1}","PEL");
	 if (Mbbb)   legendDMC->AddEntry(grdataMbbb,               "ATLAS data, 4.6 fb^{-1}","PEL");
	 if (Nbbb)   legendDMC->AddEntry(grdataNbbb,               "ATLAS data, 4.6 fb^{-1}","PEL");	 
	 legendDMC->Draw();
	 if (hname=="j_bjet_ZdPhiF20_binned") 
	 {
	     ATLAS_LABEL(0.2,0.55);
	     myText(0.33,0.55,1,"Internal");
	     myText(0.65,0.15,1,"Z+b selection");	     
	 }
	 else if (hname=="j_bjet_ZdR620_binned") 
	 {
	     ATLAS_LABEL(0.2,0.55);
	     myText(0.33,0.55,1,"Internal");
	     myText(0.65,0.25,1,"Z+b selection");	     
	 }
	 else 
	 {
	     ATLAS_LABEL(0.2,0.87);
	     myText(0.33,0.87,1,"Internal");
	     myText(0.6,0.5,1,"Z+b selection");
	 }
	 

	 c2DandRPDF->cd();

	   // rat0D = getGraphRatio(grdataMBayes,grdataStat);
	   // rat1D = getGraphRatio(grdataMBayes,grdataSyst);
	   // rat3D = getGraphRatio(grdataMBayes,grMCFM8NPsta_as_pdf_sca_np);
	 grMCFM8NPsta->SetMarkerSize(0.01);
	 grMCFM8NPsta->SetLineColor(4);
	 grMCFM8NPsta->SetLineWidth(2);
	 TGraphAsymmErrors* rat3Dstat = getGraphRatio(grdataMBayes,grMCFM8NPsta);
	 TGraphAsymmErrors* rat3Dstat_as = getGraphRatio(grdataMBayes,grMCFM8NPsta_as);
	 TGraphAsymmErrors* rat3Dstat_as_pdf = getGraphRatio(grdataMBayes,grMCFM8NPsta_as_pdf);
	 TGraphAsymmErrors* rat3Dstat_as_pdf_sca = getGraphRatio(grdataMBayes,grMCFM8NPsta_as_pdf_sca);
	 grMCFM8NPsta_as_pdf_sca_np->SetFillColor(colfMCFMall);
	 grMCFM8NPsta_as_pdf_sca_np->SetFillStyle(1001);
	 TGraphAsymmErrors* rat3Dstat_as_pdf_sca_np = getGraphRatio(grdataMBayes,grMCFM8NPsta_as_pdf_sca_np);
	   // rat4D = getGraphRatio(grdataMBayes,grMCFM23NPsta_as_pdf_sca_np);
	   // rat5D = getGraphRatio(grdataMBayes,graMCNLOMPIfb);
	   // rat6D = getGraphRatio(grdataMBayes,grSherpaNPfb);

	 hemptyRat->Draw();
	 rat0D->SetFillColor(21);
	 rat0D->Draw("E2");
	 rat1D->Draw("E2SAME");
	 //rat2D->Draw("EPCSAME");
	 //rat4D->Draw("EPCSAME");
	 rat3Dstat_as_pdf_sca_np->Draw("E2SAME");
	 rat3Dstat_as_pdf_sca->Draw("E2SAME");
	 rat3Dstat_as_pdf->Draw("E2SAME");
	 rat3Dstat_as->Draw("E2SAME");
	 rat3Dstat->Draw("E2SAME");
	 rat3Dstat->Draw("EPSAME");


	 TLegend* legendRATDpdf;
	 TLegend* legendRATSpdf; 
	 TLegend* legendRATSpdf1;

	 if (legendRATDpdf==0) {
	   //if( "j_bjet_ZdR620_binned"==hname) legendRATD=new TLegend(0.175,0.55,0.475,0.90);
	   legendRATDpdf=new TLegend(0.17,0.80,0.5,0.97);  
	 }
	 legendRATDpdf->SetShadowColor(1);
	 legendRATDpdf->SetBorderSize(0);
	 legendRATDpdf->SetFillStyle(0);
	 //legend->SetFillColor(kWhite);
	 legendRATDpdf->SetTextSize(0.06);  
	 if (legendRATSpdf==0) {
	   //legendRATS=new TLegend(0.5,0.70,0.98,0.97);   // with 4 theories 
	   legendRATSpdf=new TLegend(0.5,0.65,0.98,0.97);   // with 2 theories 
	 }
	 legendRATSpdf->SetShadowColor(1);
	 legendRATSpdf->SetBorderSize(0);
	 legendRATSpdf->SetFillStyle(0);
	 //legend->SetFillColor(kWhite);
	 legendRATSpdf->SetTextSize(0.06);  
	 if (legendRATSpdf1==0) {
	   //legendRATS=new TLegend(0.5,0.70,0.98,0.97);   // with 4 theories 
	   legendRATSpdf1=new TLegend(0.5,0.75,0.98,0.97);   // with 2 theories 
	 }
	 legendRATSpdf1->SetShadowColor(1);
	 legendRATSpdf1->SetBorderSize(0);
	 legendRATSpdf1->SetFillStyle(0);
	 //legend->SetFillColor(kWhite);
	 legendRATSpdf1->SetTextSize(0.06);  





	 legendRATDpdf->Clear();
 	 legendRATDpdf->AddEntry(rat0D,"Data stat+syst err.","F");
 	 legendRATDpdf->AddEntry(rat1D,"Data syst err.","F");
	 legendRATSpdf->Clear();
	 legendRATSpdf->AddEntry(rat3Dstat,                "MCFM(mstw8)xNP, stat","PEL");
	 legendRATSpdf->AddEntry(rat3Dstat_as,             "MCFM(mstw8)xNP, and #alpha_{s}","F");
	 legendRATSpdf->AddEntry(rat3Dstat_as_pdf,         "MCFM(mstw8)xNP, and pdf","F");
	 legendRATSpdf->AddEntry(rat3Dstat_as_pdf_sca,     "MCFM(mstw8)xNP, and scale","F");
	 legendRATSpdf->AddEntry(rat3Dstat_as_pdf_sca_np,  "MCFM(mstw8)xNP, and NP","F");
	 //legendRATSpdf->AddEntry(rat4D,"MCFM(nnpdf23)xNP, stat+syst err. ","L");
	 //legendRATS->AddEntry(rat5D,"aMC@NLO","L");
	 //legendRATS->AddEntry(rat6D,"SHERPA","L");
	 legendRATDpdf->Draw();
	 legendRATSpdf->Draw();


	 c3DandRPDF->cd();


	 hemptyRat->Draw();
	 rat0D->Draw("E2");
	 rat1D->Draw("E2SAME");
	 //rat2D->Draw("EPSAME");
	 rat2D->Draw("EPCSAME");
	 rat4D->Draw("EPCSAME");
	 rat3D->Draw("EPCSAME");



	 //rat4D->Draw("EPSAME");
	 legendRATSpdf1->Clear();
	 legendRATSpdf1->AddEntry(rat3D,"MCFM(mstw8)xNP, stat+syst err.","L");
	 legendRATSpdf1->AddEntry(rat2D,"MCFM(ct10)xNP, stat+syst err. ","L");
	 legendRATSpdf1->AddEntry(rat4D,"MCFM(nnpdf23)xNP, stat+syst err. ","L");
	 //legendRATS->AddEntry(rat5D,"aMC@NLO","L");
	 //legendRATS->AddEntry(rat6D,"SHERPA","L");
	 legendRATDpdf->Draw();
	 legendRATSpdf1->Draw();
	 
	 if (saveAsEps)
	   {
	       if (MBayes) fname = "fig"+hname+"_MCFM_w_wo_NP_dataPDF_WITHRATIO.eps";
	       else if (NBayes) fname = "fig"+hname+"_MCFM_w_wo_NP_dataPDFNic_WITHRATIO.eps"; 
	     if (setLogY) fname = "log"+fname;
	     cDandRPDF->SaveAs(fname.c_str());
	   }

	 //	 return;
	 

	 TCanvas* compareNM = new TCanvas("Comparison","Comparison",600,600);
	 hemptyRatF->Draw();
	 //rat1D->Draw("PX");
	 TGraphAsymmErrors* rat2DNM=getGraphRatio(grdataMBayes,grdataMbbb);
	 TGraphAsymmErrors* rat3DNM=getGraphRatio(grdataMBayes,grdataNbbb);
	 TGraphAsymmErrors* rat4DNM=getGraphRatio(grdataMBayes,grdataNBayes);
	 rat2DNM->Draw("PESAME");
	 rat3DNM->Draw("PESAME");
	 rat4DNM->Draw("PESAME");
// 	 legendRATD->Clear();
// 	 legendRATD->AddEntry(rat2DNM,"Mark bin by byn","PL") ;
// 	 legendRATD->AddEntry(rat4DNM,"Nicola Bayesian","PL") ;
// 	 legendRATD->AddEntry(rat3DNM,"Nicola bin by byn","PL"); 
// 	 legendRATD->Draw();
 	 legend->Clear();
 	 legend->AddEntry(rat2DNM,"An. (1) bin by byn","PL") ;
 	 legend->AddEntry(rat4DNM,"An. (2) Bayesian","PL") ;
 	 legend->AddEntry(rat3DNM,"An. (2) bin by byn","PL"); 
 	 legend->Draw();
	 if (saveAsEps)
	   {
	     fname = "fig"+hname+"_RATIO_dataChecks.eps";
	     if (setLogY) fname = "log"+fname;
	     compareNM->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
*///////////////