#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"

/// .x makeDistribPlotWratio_v3syst_trynew.C("j_bjet_Pt_binned", false, true)
/// .x makeDistribPlotWratio_v3syst_trynew.C("j_bjet_Y_binned", false, false)
/// .x makeDistribPlotWratio_v3syst_trynew.C("j_bjet_Yboost_binned", false, false)
/// .x makeDistribPlotWratio_v3syst_trynew.C("j_bjet_DY20_binned", false, true)
/// .x makeDistribPlotWratio_v3syst_trynew.C("j_bjet_ZdPhiF20_binned", true, false)
/// .x makeDistribPlotWratio_v3syst_trynew.C("j_bjet_ZdR620_binned", true, false)
/// .x makeDistribPlotWratio_v3syst_trynew.C("j_bjet_ZPt_binned", false, true)
/// .x makeDistribPlotWratio_v3syst_trynew.C("j_bjet_ZY_binned", false, false)

#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
#define showZbInPlots false
#define showLines true


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

    bool useOldaMCNLO=true;
    
  bool oldaMCNLO = useOldaMCNLO;
  bool skipMCFM = false;
  bool showWithLinSum = true;
    
  double gFSRcorr = 0.97;
  double err_gFSRcorr = 0.01;
  double err2_gFSRcorr = err_gFSRcorr*err_gFSRcorr;

  
  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 - Aug 22nd
  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_Yboost_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_bjet_DY20_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_ZY_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_ZPt_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_ZPtv3_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_dR_binned\""<<std::endl; 
    std::cout<<"\"j_2bjet_Mbb_binned\""<<std::endl;

    return;
  }
  //if (hname=="j_bjet_DY20_binned" || hname=="j_bjet_Yboost_binned" ) skipMCFM = true;
     
     TFile* fileMarkAug19comb = new TFile(  "../data/fromMarkSep2/unfolded_combined_unfold_MV175_MET70_combSF_fit_MV175_MET70_combSF.root");
     //     TFile* fileMarkAug19comb = new TFile(  "../data/fromMarkAug19/unfolded_combined_MV175_MET30_combSF.root");
     TFile* fileMarkAug19ee   = new TFile(  "../data/fromMarkAug19/unfolded_egamma_MV175_MET30_combSF.root");
     TFile* fileMarkAug19mm   = new TFile(  "../data/fromMarkAug19/unfolded_muons_MV175_MET30_combSF.root");

     TFile* fileSHPY = new TFile(  "../SherpaPythia/SherpaAndPythiaResultsNP_v3.root");
     TFile* fileaMCNLOMPI; 
     //     fileaMCNLOMPI = new TFile("../aMCatNLO/Zbb_allMergedToMPI.v3.Rivet.root");
     fileaMCNLOMPI = new TFile("../aMCatNLO/LouvainProd/MergedToMPI_goodScaling_histo_Louvain_4fZbb_MCFMscales.v3.Rivet.root");
     TFile* fileaMCNLO   = new TFile("../aMCatNLO/Zbb_all.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* 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/lecceProd/histo_oldSet4fZbb.v3.Rivet.root");
     TFile* fileaMCNLO5fZb_v3      = new TFile("../aMCatNLO/lecceProd/histo_oldSet5fZb.v3.Rivet.root");
     //TFile* fileaMCNLO5fZb_MI_v3   = new TFile("../aMCatNLO/lecceProd/histo_oldSet5fZb_MI.v3.Rivet.root");
     //     TFile* fileaMCNLO5fZb_MI_v3   = new TFile("../aMCatNLO/lecceProd/histo_oldSetNewCuts_MCFMscales_5fZb.v3.Rivet.root");
     TFile* fileaMCNLO5fZb_MI_v3   = new TFile("../aMCatNLO/lecceProd/MergedToMPI_goodScaling_histo_oldSetNewCuts_MCFMscales_5fZbMPI_1_6samples.v3.Rivet.root");
     TFile* fileaMCNLO5fZj_MI_v3   = new TFile("../aMCatNLO/LouvainProd/MergedToMPI_goodScaling_histo_Louvain_5fZj_MCFMscales_1to5.v3.Rivet.root");
     //TFile* fileaMCNLO5fZj_MI_v3   = new TFile("../aMCatNLO/LouvainProd/MergedToMPI_goodScaling_histo_Louvain_5fZj_MCFMscales_1to10.v3.Rivet.root");
     //TFile* fileaMCNLO5fZj_MI_v3   = new TFile("../aMCatNLO/lecceProd/MergedToMPI_histo_oldSetNewCuts_MCFMscales_5fZj_1_2samples.v3.Rivet.root");
     //TFile* fileaMCNLO5fZb_MI_v3   = new TFile("../aMCatNLO/lecceProd/histo_oldSetNewCuts_Aug22_5fZb_wMPI100k_buggyscales.v3.Rivet.root");
     TFile* fileaMCNLO5fZb_noFR_v3 = new TFile("../aMCatNLO/lecceProd/histo_oldSet5fZb_noFR.v3.Rivet.root");
     TFile* fileaMCNLO5fZj_v3      = new TFile("../aMCatNLO/lecceProd/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;

     TFile* fileNikoTotSyst = NULL;

     if (hname.substr(0,3)=="j_b") 
     {
	 fileMCFM      = new TFile("../MCFMall/ZbMCFMgraphCT10_aug13.root");
	 fileMCFMmstw  = new TFile("../MCFMall/ZbMCFMgraphMSTW2008_aug13.root");
	 fileMCFMnnpdf = new TFile("../MCFMall/ZbMCFMgraphNNPDF23_aug13.root");
     }
     else if (hname.substr(0,3)=="j_2")
     {
	 fileMCFM      = new TFile("../MCFMall/ZbbMCFMgraphCT10_aug13.root");
	 fileMCFMmstw  = new TFile("../MCFMall/ZbbMCFMgraphMSTW2008_aug13.root");
	 fileMCFMnnpdf = new TFile("../MCFMall/ZbbMCFMgraphNNPDF_aug13.root");
	 fileAlpgen = 0;
     }




     // 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 *graMCNLO5fZj_MI; 
     TGraphAsymmErrors *graMCNLO5fZb_noFR; 
     TGraphAsymmErrors *grSherpafb; // in fb 
     TGraphAsymmErrors *grPythiafb; // in fb 
     TGraphAsymmErrors *graMCNLOfb; 
     TGraphAsymmErrors *graMCNLO4ffb; 
     TGraphAsymmErrors *graMCNLO4flsfb; 
     TGraphAsymmErrors *graMCNLO5fZbfb; 
     TGraphAsymmErrors *graMCNLO5fZjfb; 
     TGraphAsymmErrors *graMCNLO5fZj_MIfb; 
     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 MCFMNP_EYl_totquadrsum[100];
     Double_t MCFMNP_EYh_totquadrsum[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 MCFM8NP_EYl_totquadrsum[100];
     Double_t MCFM8NP_EYh_totquadrsum[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];
     Double_t MCFM23NP_EYl_totquadrsum[100];
     Double_t MCFM23NP_EYh_totquadrsum[100];


     Double_t MCFM_errNP_H[100];
     Double_t MCFM_errNP_L[100];
     Double_t MCFM8_errNP_H[100];
     Double_t MCFM8_errNP_L[100];
     Double_t MCFM23_errNP_H[100];
     Double_t MCFM23_errNP_L[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 > 0) std::cout<<" fdataname will be "<<fdataname<<std::endl;
     

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

     TH2F* hempty; 
     TH2F* hemptyRat; 
     TH2F* hemptyRatPDF; 
     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;
     std::string hnameMarkAug19;
     int offsetMarkAug19=0;
     //     TFile* fileMarkAug19comb = new TFile(  "../data/fromMarkAug19/unfolded_combined_MV175_MET30_combSF.root.root");
     //     TFile* fileMarkAug19ee = new TFile(  "../data/fromMarkAug19/unfolded_egamma_MV175_MET30_combSF.root.root");
     //     TFile* fileMarkAug19mm = new TFile(  "../data/fromMarkAug19/unfolded_muons_MV175_MET30_combSF.root.root");
     if (hname=="j_bjet_Pt_binned")
       {
         fileNikoTotSyst = new TFile(  "../data/fromNikoTotSystSept2/CombinedSystematics_JetPt.root");
	 hnameMarkAug19=  "bjetpt_unfolded_bayes_h";	   
	 offsetMarkAug19=2;	 
	 fileAlpgen= new TFile("../alpgen/newV3/AlpgenUnfoldingSummary_CombinedChannel_JetPt_Central.root");
	 xTitleStr = "(b-jet selection) b-jet p_{T} [GeV]";
	 xTitleStr = "b-jet p_{T} [GeV]";
	 yTitleStr = "d#sigma/dp_{T} [fb/GeV]";
	 hnameMCFM = "Pt_gae_";
	 hnameaMCNLO = "bjet_pT";
	 hnameaMCNLOmpi = hname;

	 double xmin=0.;
	 double xmax=600.;
	 double ymin=0.;
	 double ymax=250.;
	 double yminRat=0.3;
	 double ymaxRat=1.7;
	 double yminRatPDF=yminRat;
	 double ymaxRatPDF=ymaxRat;
	 if (setLogX) 
	 {
	     xmin=17.;
	     ymin=-5.;
	     xmax=600.;
	 }
	 if (setLogY)
	 {
	     ymin=0.05;
	     ymax=1000.;
	 }	 
	 hempty     = new TH2F("empty","empty"        ,10,xmin,xmax,10,ymin   ,ymax   ); 
	 hemptyRat  = new TH2F("emptyRat","emptyRat"  ,10,xmin,xmax,10,yminRat,ymaxRat); 
	 hemptyRatPDF  = new TH2F("emptyRatPDF","emptyRatPDF"  ,10,xmin,xmax,10,yminRatPDF,ymaxRatPDF); 
	 hemptyRatF = new TH2F("emptyRatF","emptyRatF",10,xmin,xmax,10,ratfmin,ratfmax); 

	 //grEmpty->GetXaxis()->SetRangeUser(0.,600.);
       }
     else if ("j_bjet_Y_binned"==hname)
       {
         fileNikoTotSyst = new TFile(  "../data/fromNikoTotSystSept2/CombinedSystematics_JetY.root");
	 hnameMarkAug19=  "bjetrap_unfolded_bayes_h";	   
	 fileAlpgen= new TFile("../alpgen/newV3/AlpgenUnfoldingSummary_CombinedChannel_JetY_Central.root");
	 xTitleStr = "(b-jet selection) b-jet |Y|";
	 xTitleStr = "b-jet |y|";
	 yTitleStr = "d#sigma/d|y| [fb/y-unit]";

	 double xmin=0.;
	 double xmax=2.5;
	 double ymin=400.;
	 double ymax=3750.;
	 double yminRat=0.45;
	 double ymaxRat=1.55;
	 double yminRatPDF=yminRat;
	 double ymaxRatPDF=ymaxRat;
	 if (setLogX) 
	 {
	     xmin=17.;
	     ymin=-5.;
	     xmax=600.;
	 }
	 if (setLogY)
	 {
	     ymin=0.05;
	     ymax=1000.;
	 }
	 hempty     = new TH2F("empty","empty"        ,10,xmin,xmax,10,ymin   ,ymax   ); 
	 hemptyRat  = new TH2F("emptyRat","emptyRat"  ,10,xmin,xmax,10,yminRat,ymaxRat); 
	 hemptyRatPDF  = new TH2F("emptyRatPDF","emptyRatPDF"  ,10,xmin,xmax,10,yminRatPDF,ymaxRatPDF); 
	 hemptyRatF = new TH2F("emptyRatF","emptyRatF",10,xmin,xmax,10,ratfmin,ratfmax); 

	 hnameMCFM = "Y_gae_";
	 hnameaMCNLO = "bjet_y";
	 hnameaMCNLOmpi = hname;
       }
     else if ("j_bjet_Yboost_binned"==hname)
       {
         fileNikoTotSyst = new TFile(  "../data/fromNikoTotSystSept2/CombinedSystematics_Ystar.root");
	 hnameMarkAug19=  "bjetZyboost_unfolded_bayes_h";	   
	 fileAlpgen= new TFile("../alpgen/newV3/AlpgenUnfoldingSummary_CombinedChannel_Ystar_Central.root");
	 xTitleStr = "(b-jet selection) b-jet Y_{boost}";
	 xTitleStr = "y_{boost}";
	 yTitleStr = "d#sigma/dy_{boost} [fb/y-unit]";

	 double xmin=0.;
	 double xmax=2.6;
	 double ymin=-100.;
	 double ymax=5000.;
	 double yminRat=0.45;
	 double ymaxRat=1.55;
	 double yminRatPDF=yminRat;
	 double ymaxRatPDF=ymaxRat;
	 if (setLogX) 
	 {
	     xmin=17.;
	     ymin=-5.;
	     xmax=600.;
	 }
	 if (setLogY)
	 {
	     ymin=0.05;
	     ymax=1000.;
	 }
	 hempty     = new TH2F("empty","empty"        ,10,xmin,xmax,10,ymin   ,ymax   ); 
	 hemptyRat  = new TH2F("emptyRat","emptyRat"  ,10,xmin,xmax,10,yminRat,ymaxRat); 
	 hemptyRatPDF  = new TH2F("emptyRatPDF","emptyRatPDF"  ,10,xmin,xmax,10,yminRatPDF,ymaxRatPDF); 
	 hemptyRatF = new TH2F("emptyRatF","emptyRatF",10,xmin,xmax,10,ratfmin,ratfmax); 

	 hnameMCFM = "Yboost_gae_";
	 hnameaMCNLO = hname;
	 hnameaMCNLOmpi = hname;
       }
     else if ("j_bjet_ZY_binned"==hname)
       {
         fileNikoTotSyst = new TFile(  "../data/fromNikoTotSystSept2/CombinedSystematics_ZY.root");
	 hnameMarkAug19=  "Zrap_unfolded_bayes_h";	   
	 fileAlpgen= new TFile("../alpgen/newV3/AlpgenUnfoldingSummary_CombinedChannel_ZY_Central.root");
	 xTitleStr = "(Zb(b) selection Z |Y|";
	 xTitleStr = "Z |y|";
	 yTitleStr = "d#sigma/d|y| [fb/y-unit]";

	 double xmin=0.;
	 double xmax=2.6;
	 double ymin=-100.;
	 double ymax=4500.;
	 double yminRat=0.45;
	 double ymaxRat=1.55;
	 double yminRatPDF=yminRat;
	 double ymaxRatPDF=ymaxRat;
	 if (setLogX) 
	 {
	     xmin=17.;
	     ymin=-5.;
	     xmax=600.;
	 }
	 if (setLogY)
	 {
	     ymin=0.05;
	     ymax=1000.;
	 }
	 hempty     = new TH2F("empty","empty"        ,10,xmin,xmax,10,ymin   ,ymax   ); 
	 hemptyRat  = new TH2F("emptyRat","emptyRat"  ,10,xmin,xmax,10,yminRat,ymaxRat); 
	 hemptyRatPDF  = new TH2F("emptyRatPDF","emptyRatPDF"  ,10,xmin,xmax,10,yminRatPDF,ymaxRatPDF); 
	 hemptyRatF = new TH2F("emptyRatF","emptyRatF",10,xmin,xmax,10,ratfmin,ratfmax); 

	 hnameMCFM = "ZY_gae_";
	 hnameaMCNLO = "UFFA";
	 hnameaMCNLOmpi = hname;
       }
     else if ("j_bjet_ZPt_binned"==hname)
       {
	 hnameMarkAug19=  "Zpt_unfolded_bayes_h";	   
         fileNikoTotSyst = new TFile(  "../data/fromNikoTotSystSept2/CombinedSystematics_ZPt.root");
	 fileAlpgen= new TFile("../alpgen/newV3/AlpgenUnfoldingSummary_CombinedChannel_ZPt_Central.root");
	 xTitleStr = "Z p_{T} [GeV]";
	 yTitleStr = "d#sigma/dp_{T} [fb/GeV]";

	 double xmin=-20.;
	 double xmax=550.;
	 double ymin=-5.;
	 double ymax=125.;
	 double yminRat=0.25;
	 double ymaxRat=1.95;
	 double yminRatPDF=yminRat;
	 double ymaxRatPDF=ymaxRat;
	 if (setLogX) 
	 {
	     xmin=3.;
	     xmax=1000.;
	     ymin=-5.;
	 }
	 if (setLogY)
	 {
	     ymin=0.01;
	     ymax=400.;
	 }	 
	 hempty     = new TH2F("empty","empty"        ,10,xmin,xmax,10,ymin   ,ymax   ); 
	 hemptyRat  = new TH2F("emptyRat","emptyRat"  ,10,xmin,xmax,10,yminRat,ymaxRat); 
	 hemptyRatPDF  = new TH2F("emptyRatPDF","emptyRatPDF"  ,10,xmin,xmax,10,yminRatPDF,ymaxRatPDF); 
	 hemptyRatF = new TH2F("emptyRatF","emptyRatF",10,xmin,xmax,10,ratfmin,ratfmax); 

	 hnameMCFM = "ZPt_gae_";
	 hnameaMCNLO = "Z_pT";
	 hnameaMCNLOmpi = hname;
       }
     else if ("j_bjet_ZdR620_binned"==hname)
       {
         fileNikoTotSyst = new TFile(  "../data/fromNikoTotSystSept2/CombinedSystematics_DR.root");
	 hnameMarkAug19=  "bjetZdR_unfolded_bayes_h";	   
	 fileAlpgen= new TFile("../alpgen/newV3/AlpgenUnfoldingSummary_CombinedChannel_DR_Central.root");
	 xTitleStr = "(b-jet selection) #DeltaR(Z,b-jet)";
	 xTitleStr = "#DeltaR(Z,b-jet)";
	 yTitleStr = "d#sigma/dR [fb/R-unit]";

	 double xmin=0.;
	 double xmax=6.2;
	 double ymin=-100.;
	 double ymax=5000.;
	 double yminRat=0.1;
	 double ymaxRat=1.7;
	 double yminRatPDF=yminRat;
	 double ymaxRatPDF=ymaxRat;
	 if (setLogX) 
	 {
	     xmin=3.;
	     xmax=1000.;
	     ymin=-5.;
	 }
	 if (setLogY)
	 {
	     ymin=2.;
	     ymax=7000.;
	 }	 
	 hempty     = new TH2F("empty","empty"        ,10,xmin,xmax,10,ymin   ,ymax   ); 
	 hemptyRat  = new TH2F("emptyRat","emptyRat"  ,10,xmin,xmax,10,yminRat,ymaxRat); 
	 hemptyRatPDF  = new TH2F("emptyRatPDF","emptyRatPDF"  ,10,xmin,xmax,10,yminRatPDF,ymaxRatPDF); 
	 hemptyRatF = new TH2F("emptyRatF","emptyRatF",10,xmin,xmax,10,ratfmin,ratfmax); 








// 	 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";
	 hnameaMCNLOmpi = hname;	
       }
     else if ("j_bjet_ZdPhiF20_binned"==hname)
       {
         fileNikoTotSyst = new TFile(  "../data/fromNikoTotSystSept2/CombinedSystematics_DPhi.root");
	 hnameMarkAug19=  "bjetZdphi_unfolded_bayes_h";	   	 
	 fileAlpgen= new TFile("../alpgen/newV3/AlpgenUnfoldingSummary_CombinedChannel_DPhi_Central.root");
	 xTitleStr = "(b-jet selection)  #Delta#phi(Z,b-jet)";
	 xTitleStr = "#Delta#phi(Z,b-jet)";
	 yTitleStr = "d#sigma/d#phi [fb/rad]";

	 double xmin=0.;
	 double xmax=3.3;
	 double ymin=-100.;
	 double ymax=10000.;
	 double yminRat=0.1;
	 double ymaxRat=1.7;
	 double yminRatPDF=yminRat;
	 double ymaxRatPDF=ymaxRat;
	 if (setLogX) 
	 {
	     xmin=3.;
	     xmax=1000.;
	     ymin=-5.;
	 }
	 if (setLogY)
	 {
	     ymin=50.;
	     ymax=15000.;
	 }	 
	 hempty     = new TH2F("empty","empty"        ,10,xmin,xmax,10,ymin   ,ymax   ); 
	 hemptyRat  = new TH2F("emptyRat","emptyRat"  ,10,xmin,xmax,10,yminRat,ymaxRat); 
	 hemptyRatPDF  = new TH2F("emptyRatPDF","emptyRatPDF"  ,10,xmin,xmax,10,yminRatPDF,ymaxRatPDF); 
	 hemptyRatF = new TH2F("emptyRatF","emptyRatF",10,xmin,xmax,10,ratfmin,ratfmax); 


	 hnameMCFM = "ZdPhi_gae_";
	 hnameaMCNLO = "Zj_dphi";
	 hnameaMCNLOmpi = hname;
       }
     else if ("j_bjet_DY20_binned"==hname)
       {
         fileNikoTotSyst = new TFile(  "../data/fromNikoTotSystSept2/CombinedSystematics_DY.root");
	 hnameMarkAug19=  "bjetZdy_unfolded_bayes_h";	   
	 fileAlpgen= new TFile("../alpgen/newV3/AlpgenUnfoldingSummary_CombinedChannel_DY_Central.root");
	 xTitleStr = "(b-jet selection) |#Delta y(Z,b-jet)|";
	 xTitleStr = "|#Delta y(Z,b-jet)|";
	 yTitleStr = "d#sigma/d|#Delta y| [fb/Y-unit]";
	 double xmin=0.;
	 double xmax=2.5;
	 double ymin=400.;
	 double ymax=4000.;
	 double yminRat=0.35;
	 double ymaxRat=1.55;
	 double yminRatPDF=yminRat;
	 double ymaxRatPDF=ymaxRat;
	 if (setLogX) 
	 {
	     xmin=0.06;
	     xmax=7.;
	     ymin=-100.;
	 }
	 if (setLogY)
	 {
	     ymin=0.05;
	     ymax=1000.;
	 }
	 hempty     = new TH2F("empty","empty"        ,10,xmin,xmax,10,ymin   ,ymax   ); 
	 hemptyRat  = new TH2F("emptyRat","emptyRat"  ,10,xmin,xmax,10,yminRat,ymaxRat); 
	 hemptyRatPDF  = new TH2F("emptyRatPDF","emptyRatPDF"  ,10,xmin,xmax,10,yminRatPDF,ymaxRatPDF); 
	 hemptyRatF = new TH2F("emptyRatF","emptyRatF",10,xmin,xmax,10,ratfmin,ratfmax); 

	 hnameMCFM = "ZdY_gae_";
	 hnameaMCNLO = "UFFA";
	 hnameaMCNLOmpi = hname;
	 //	 std::cout<<" siamo entrati qui"<<std::endl;
	 
       }
     else if ("j_2bjet_ZPt_binned"==hname || "j_2bjet_ZPtv3_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";
	 hnameaMCNLOmpi = hname;
       }
     else if ("j_2bjet_ZY_binned"==hname)
       {
	 hnameMarkAug19=  "XXX";	   
	 xTitleStr = "(Zbb selection) Z |Y|";
	 yTitleStr = "d#sigma/d|Y| [fb/Y-unit]";
	 if (setLogY) hempty = new TH2F("empty","empty",10,0.0,3.,10,100.,500.0); // axis range hpx->SetStats(kFALSE);	// no statistics
	 else hempty = new TH2F("empty","empty",10,0.0,3.,10,.0,500.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 = hname;
	 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";
	 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";
	 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 ) legendDMC=new TLegend(0.175,0.65,0.475,0.90);
       else if(  "j_bjet_ZdR620_binned"==hname   ) legendDMC=new TLegend(0.2,0.03,0.5,0.28);
       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.65,0.8,1.,0.97);
	 else 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 
	 if( "j_bjet_ZdR620_binned"==hname) 	legendRATS=new TLegend(0.17,0.75,0.65,0.97);   
	 else 	legendRATS=new TLegend(0.5,0.75,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 
     
     TGraphAsymmErrors* grMarkAug19comb = NULL;
     TGraphAsymmErrors* grMarkAug19combWithSystErrs;

     TH1F *htemp = (TH1F*)fileMarkAug19comb->Get(hnameMarkAug19.c_str());
     if (htemp!=NULL)  grMarkAug19comb = new TGraphAsymmErrors(htemp);
     
     //if (fileNikoTotSyst!=NULL) htemp = (TH1F*)fileNikoTotSyst->Get("Tot_up_Temp");
     
     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());
     graMCNLO5fZj_MI   = (TGraphAsymmErrors*)fileaMCNLO5fZj_MI_v3  ->Get(hname.c_str());
     graMCNLO5fZb_noFR = (TGraphAsymmErrors*)fileaMCNLO5fZb_noFR_v3->Get(hname.c_str());
     graMCNLOMPI       = (TGraphAsymmErrors*)fileaMCNLOMPI         ->Get(hname.c_str());
     graMCNLO          = (TGraphAsymmErrors*)fileaMCNLO            ->Get(hname.c_str());
     //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;        graMCNLOMPIX      = graMCNLOMPI->GetX();     
     Double_t* graMCNLOMPIY;        graMCNLOMPIY      = graMCNLOMPI->GetY();
     Double_t* graMCNLOMPIEyLow;    graMCNLOMPIEyLow  = graMCNLOMPI->GetEYlow();
     Double_t* graMCNLOMPIEyHigh;   graMCNLOMPIEyHigh = graMCNLOMPI->GetEYhigh();
     Double_t* graMCNLOX;        graMCNLOX      = graMCNLO->GetX();
     Double_t* graMCNLOY;        graMCNLOY      = graMCNLO->GetY();
     Double_t* graMCNLOEyLow;    graMCNLOEyLow  = graMCNLO->GetEYlow();
     Double_t* graMCNLOEyHigh;   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* graMCNLO5fZj_MIY;        graMCNLO5fZj_MIY      = graMCNLO5fZj_MI->GetY();
     Double_t* graMCNLO5fZj_MIEyLow;    graMCNLO5fZj_MIEyLow  = graMCNLO5fZj_MI->GetEYlow();
     Double_t* graMCNLO5fZj_MIEyHigh;   graMCNLO5fZj_MIEyHigh = graMCNLO5fZj_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  graMCNLO5fZj_MIYfb[100]; 
     Double_t  graMCNLO5fZj_MIEyLowfb[100];
     Double_t  graMCNLO5fZj_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;
      int Npoints = grSherpa->GetN();
      if (hname.substr(0,6)=="j_bjet") {
	/*
	readVectorsFromFile(fdataname,nrows,X1,X2,EX2,X3,X4,EX4,sUP,sDW,m_debug);
	  if (nrows==0) 
	  {
	      for (int j=0; j<grMarkAug19comb->GetN(); ++j)
	      {
		  dataMBayes[j]  = grMarkAug19comb->GetY()[j+offsetMarkAug19];
		  edataMBayes[j] = grMarkAug19comb->GetEYlow()[j+offsetMarkAug19];
		  dataMbbb[j]    = dataMBayes[j];

		  dataNbbb[j]   = dataMBayes[j];
		  dataNBayes[j] = dataMBayes[j];
		  edataNBayes[j]= edataMBayes[j];

		  // fake syst error = 10% 
		  edataMBayesSup[j]= 0.1*dataMBayes[j];
		  edataMBayesSdw[j]= 0.1*dataMBayes[j];	   
		  edataNBayesSup[j] = 0.1*dataNBayes[j];
		  edataNBayesSdw[j] = 0.1*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;		  
	      }
	  }
	  else 
	  {
	      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];
		  dataMBayes[j] = grMarkAug19comb->GetY()[j+offsetMarkAug19];
		  //edataMBayes[j] = EX4[j];
		  edataMBayes[j] = grMarkAug19comb->GetEYlow()[j+offsetMarkAug19];

		  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;
	      }
	  }
	  */
	TH1F* htempSysUP=NULL;
	TH1F* htempSysSW=NULL;
	if (fileNikoTotSyst!=NULL) {
	  htempSysUP = (TH1F*)fileNikoTotSyst->Get("Tot_up_Temp");
	  htempSysDW = (TH1F*)fileNikoTotSyst->Get("Tot_down_Temp");
	  TCanvas * tempSyst = new TCanvas("syst","syst");
	  htempSysUP->Draw("C");
	  htempSysDW->Draw("CSAME");
	}
	else {
	  std::cout<<"There's a problem: data Total syst not available"<<std::endl;
	  return;
	}
	  for (int j=0; j<Npoints; ++j)
	    {
	      dataNbbb[j]   = X1[j];
	      dataNBayes[j] = X2[j];
	      edataNBayes[j]= EX2[j];
	      
	      dataMBayes[j] = grMarkAug19comb->GetY()[j+offsetMarkAug19];
	      dataMbbb[j]   = dataMBayes[j];
	      edataMBayes[j] = grMarkAug19comb->GetEYlow()[j+offsetMarkAug19];
	      edataMbbb[j]   = edataMBayes[j];
	      
	      //
	      edataMBayesSup[j]=  dataMBayes[j] * htempSysUP->GetBinContent(j+1);
	      edataMBayesSdw[j]= -dataMBayes[j] * htempSysDW->GetBinContent(j+1);
	      //

	      edataNBayesSup[j] = edataMBayesSup[j]; 
	      edataNBayesSdw[j] = edataMBayesSdw[j];
	      // edataNbbbSup[j]   = edataMBayesSup[j]; 
	      // edataNbbbSdw[j]   = edataMBayesSdw[j]; 
	      // edataMbbbSup[j]   = edataMBayesSup[j]; 
	      // edataMbbbSdw[j]   = edataMBayesSdw[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 == "<<dataMBayes[j]<<" stat, Up/Dw Syst "<<edataMBayes[j]<<"  "<<edataMBayesSup[j]<<"/"<<edataMBayesSdw[j]
<<" -------- Up/Dw Tot "<<"  "<<edataMBayesSSup[j]<<"/"<<edataMBayesSSdw[j] <<std::endl;
	    }
      }
      //return;
      if (m_debug>1) std::cout<<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% data read"<<std::endl;


     
     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]  = 0.5*sqrt(grNPSherpaEYl[k]*grNPSherpaEYl[k]+grNPPythiaEYl[k]*grNPPythiaEYl[k]); // multiplicative term
	 grNPEYhstat[k]  = 0.5*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]       = gFSRcorr* MCFMNP_Y[k];  // after FSR 
	 MCFM8NP_Y[k]      = gFSRcorr* MCFM8NP_Y[k]; // after FSR 
	 MCFM23NP_Y[k]     = gFSRcorr* MCFM23NP_Y[k];// after FSR 
	 
	 double corr = gFSRcorr * grNPY[k];
	 /////////// temporary 
	 /////////// corr = 1.;
	 ///////////


	 double staErrL = corr * grMCFMsta->GetEYlow()[k];
	 double staErrH = corr * grMCFMsta->GetEYhigh()[k];
	 double asErrL  = corr * grMCFMas ->GetEYlow()[k];
	 double asErrH  = corr * grMCFMas ->GetEYhigh()[k];
	 double pdfErrL = corr * grMCFMpdf->GetEYlow()[k];
	 double pdfErrH = corr * grMCFMpdf->GetEYhigh()[k];
	 double scaErrL = corr * grMCFMsca->GetEYlow()[k];
	 double scaErrH = corr * 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]*err2_gFSRcorr); // multiplicative errors from NP and 1% for FSR
// 	 npErrH  = sqrt(npErrH*npErrH + MCFMNP_Y[k]*MCFMNP_Y[k]*err2_gFSRcorr); // multiplicative errors from NP and 1% for FSR
	 double npErrL = MCFMNP_Y[k] * sqrt(  grNPEYlstatAndSys[k]*grNPEYlstatAndSys[k]/(grNPY[k]*grNPY[k]) + 
					      err2_gFSRcorr/(gFSRcorr*gFSRcorr) );
	 double npErrH = MCFMNP_Y[k] * sqrt(  grNPEYhstatAndSys[k]*grNPEYhstatAndSys[k]/(grNPY[k]*grNPY[k]) + 
					      err2_gFSRcorr/(gFSRcorr*gFSRcorr) );
	 MCFM_errNP_L[k]=npErrL;
	 MCFM_errNP_H[k]=npErrH;

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

         MCFMNP_EYlstat[k] = staErrL;
	 if (showWithLinSum) 
	 {
	     MCFMNP_EYlstat_as[k]            = MCFMNP_EYlstat[k]           + asErrL;
	     MCFMNP_EYlstat_as_pdf[k]        = MCFMNP_EYlstat_as[k]        + pdfErrL;
	     MCFMNP_EYlstat_as_pdf_sca[k]    = MCFMNP_EYlstat_as_pdf[k]    + scaErrL;
	     MCFMNP_EYlstat_as_pdf_sca_np[k] = MCFMNP_EYlstat_as_pdf_sca[k]+ npErrL;
	     MCFMNP_EYl_totquadrsum[k]       = sqrt(MCFMNP_EYlstat[k]*MCFMNP_EYlstat[k]                +  
					       	    asErrL*asErrL                                      +
					       	    pdfErrL*pdfErrL                                    +
					       	    scaErrL*scaErrL                                    +
						    npErrL*npErrL); 
	 }
	 else 
	 {
	     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_EYl_totquadrsum[k]       = MCFMNP_EYlstat_as_pdf_sca_np[k] ;
	 }
	 
         MCFMNP_EYhstat[k] = staErrH;
	 if (showWithLinSum) 
	 {
	     MCFMNP_EYhstat_as[k]            = MCFMNP_EYhstat[k]           + asErrH;
	     MCFMNP_EYhstat_as_pdf[k]        = MCFMNP_EYhstat_as[k]        + pdfErrH;
	     MCFMNP_EYhstat_as_pdf_sca[k]    = MCFMNP_EYhstat_as_pdf[k]    + scaErrH;
	     MCFMNP_EYhstat_as_pdf_sca_np[k] = MCFMNP_EYhstat_as_pdf_sca[k]+ npErrH;
	     MCFMNP_EYh_totquadrsum[k]       = sqrt(MCFMNP_EYhstat[k]*MCFMNP_EYhstat[k]                +  
					       	    asErrH *asErrH                                     +
					       	    pdfErrH*pdfErrH                                    +
					       	    scaErrH*scaErrH                                    +
						    npErrH *npErrH ); 
	 }
	 else 
	 {
	     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); 
	     MCFMNP_EYh_totquadrsum[k]       = MCFMNP_EYhstat_as_pdf_sca_np[k] ;
	 }
	 
	 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;
	 }

	 // MSTW2008
	 staErrL = corr * grMCFM8sta->GetEYlow()[k];
	 staErrH = corr * grMCFM8sta->GetEYhigh()[k];
	 asErrL  = corr * grMCFM8as ->GetEYlow()[k];
	 asErrH  = corr * grMCFM8as ->GetEYhigh()[k];
	 pdfErrL = corr * grMCFM8pdf->GetEYlow()[k];
	 pdfErrH = corr * grMCFM8pdf->GetEYhigh()[k];
	 scaErrL = corr * grMCFM8sca->GetEYlow()[k];
	 scaErrH = corr * grMCFM8sca->GetEYhigh()[k];
	 npErrL = MCFM8NP_Y[k] * sqrt(  grNPEYlstatAndSys[k]*grNPEYlstatAndSys[k]/(grNPY[k]*grNPY[k]) + 
					err2_gFSRcorr/(gFSRcorr*gFSRcorr) );
	 npErrH = MCFM8NP_Y[k] * sqrt(  grNPEYhstatAndSys[k]*grNPEYhstatAndSys[k]/(grNPY[k]*grNPY[k]) + 
					err2_gFSRcorr/(gFSRcorr*gFSRcorr) );
	 MCFM8_errNP_L[k]=npErrL;
	 MCFM8_errNP_H[k]=npErrH;


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

	 MCFM8NP_EYlstat[k] = staErrL;
	 if (showWithLinSum)
	 {
	     MCFM8NP_EYlstat_as[k]            = MCFM8NP_EYlstat[k]           + asErrH;
	     MCFM8NP_EYlstat_as_pdf[k]        = MCFM8NP_EYlstat_as[k]        + pdfErrH;
	     MCFM8NP_EYlstat_as_pdf_sca[k]    = MCFM8NP_EYlstat_as_pdf[k]    + scaErrH;
	     MCFM8NP_EYlstat_as_pdf_sca_np[k] = MCFM8NP_EYlstat_as_pdf_sca[k]+ npErrH;
	     MCFM8NP_EYl_totquadrsum[k]       = sqrt(MCFM8NP_EYlstat[k]*MCFM8NP_EYlstat[k]             +  
					       	    asErrL * asErrL                                    +
					       	    pdfErrL*pdfErrL                                    +
					       	    scaErrL*scaErrL                                    +
						    npErrL * npErrL ); 
	 }
	 else 
	 {
	     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_EYl_totquadrsum[k]       = MCFM8NP_EYlstat_as_pdf_sca_np[k];
	 }
         MCFM8NP_EYhstat[k] = staErrH;
	 if (showWithLinSum)
	 {
	     MCFM8NP_EYhstat_as[k]            = MCFM8NP_EYhstat[k]           + asErrH;
	     MCFM8NP_EYhstat_as_pdf[k]        = MCFM8NP_EYhstat_as[k]        + pdfErrH;
	     MCFM8NP_EYhstat_as_pdf_sca[k]    = MCFM8NP_EYhstat_as_pdf[k]    + scaErrH;
	     MCFM8NP_EYhstat_as_pdf_sca_np[k] = MCFM8NP_EYhstat_as_pdf_sca[k]+ npErrH;
	     MCFM8NP_EYh_totquadrsum[k]       = sqrt(MCFM8NP_EYhstat[k]*MCFM8NP_EYhstat[k]             +  
					       	    asErrH * asErrH                                    +
					       	    pdfErrH*pdfErrH                                    +
					       	    scaErrH*scaErrH                                    +
						    npErrH * npErrH ); 
	 }
	 else 
	 {
	     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); 
	     MCFM8NP_EYh_totquadrsum[k]       = MCFM8NP_EYhstat_as_pdf_sca_np[k];
	 }
	 
	 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;
	 }


	 // NNPDF23
	 staErrL = corr * grMCFM23sta->GetEYlow()[k];
	 staErrH = corr * grMCFM23sta->GetEYhigh()[k];
	 asErrL  = corr * grMCFM23as ->GetEYlow()[k];
	 asErrH  = corr * grMCFM23as ->GetEYhigh()[k];
	 pdfErrL = corr * grMCFM23pdf->GetEYlow()[k];
	 pdfErrH = corr * grMCFM23pdf->GetEYhigh()[k];
	 scaErrL = corr * grMCFM23sca->GetEYlow()[k];
	 scaErrH = corr * grMCFM23sca->GetEYhigh()[k];
	 npErrL = MCFM23NP_Y[k] * sqrt(  grNPEYlstatAndSys[k]*grNPEYlstatAndSys[k]/(grNPY[k]*grNPY[k]) + 
					 err2_gFSRcorr/(gFSRcorr*gFSRcorr) );
	 npErrH = MCFM23NP_Y[k] * sqrt(  grNPEYhstatAndSys[k]*grNPEYhstatAndSys[k]/(grNPY[k]*grNPY[k]) + 
					 err2_gFSRcorr/(gFSRcorr*gFSRcorr) );
	 MCFM23_errNP_L[k]=npErrL;
	 MCFM23_errNP_H[k]=npErrH;

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

	 MCFM23NP_EYlstat[k] = staErrL;
	 if (showWithLinSum)
	 {
	     MCFM23NP_EYlstat_as[k]            = MCFM23NP_EYlstat[k]           + asErrH;
	     MCFM23NP_EYlstat_as_pdf[k]        = MCFM23NP_EYlstat_as[k]        + pdfErrH;
	     MCFM23NP_EYlstat_as_pdf_sca[k]    = MCFM23NP_EYlstat_as_pdf[k]    + scaErrH;
	     MCFM23NP_EYlstat_as_pdf_sca_np[k] = MCFM23NP_EYlstat_as_pdf_sca[k]+ npErrH;
	     MCFM23NP_EYl_totquadrsum[k]       = sqrt(MCFM23NP_EYlstat[k]*MCFM23NP_EYlstat[k]          +  
					       	    asErrL * asErrL                                    +
					       	    pdfErrL*pdfErrL                                    +
					       	    scaErrL*scaErrL                                    +
						    npErrL * npErrL ); 
	 }
	 else 
	 {
	     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_EYl_totquadrsum[k]       = MCFM23NP_EYlstat_as_pdf_sca_np[k];
	 }
         MCFM23NP_EYhstat[k] = staErrH;
	 if (showWithLinSum)
	 {
	     MCFM23NP_EYhstat_as[k]            = MCFM23NP_EYhstat[k]           + asErrH;
	     MCFM23NP_EYhstat_as_pdf[k]        = MCFM23NP_EYhstat_as[k]        + pdfErrH;
	     MCFM23NP_EYhstat_as_pdf_sca[k]    = MCFM23NP_EYhstat_as_pdf[k]    + scaErrH;
	     MCFM23NP_EYhstat_as_pdf_sca_np[k] = MCFM23NP_EYhstat_as_pdf_sca[k]+ npErrH;
	     MCFM23NP_EYh_totquadrsum[k]       = sqrt(MCFM23NP_EYhstat[k]*MCFM23NP_EYhstat[k]          +  
					       	    asErrH * asErrH                                    +
					       	    pdfErrH*pdfErrH                                    +
					       	    scaErrH*scaErrH                                    +
						    npErrH * npErrH ); 
	 }
	 else 
	 {
	     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); 
	     MCFM23NP_EYh_totquadrsum[k]       = MCFM23NP_EYhstat_as_pdf_sca_np[k];
	 }

	 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];

	 graMCNLO5fZj_MIYfb[k]      = 1000.*graMCNLO5fZj_MIY[k];
	 graMCNLO5fZj_MIEyLowfb[k]  = 1000.*graMCNLO5fZj_MIEyLow[k];
	 graMCNLO5fZj_MIEyHighfb[k] = 1000.*graMCNLO5fZj_MIEyHigh[k];

	 graMCNLO5fZb_noFRYfb[k]      = 1000.*graMCNLO5fZb_noFRY[k];
	 graMCNLO5fZb_noFREyLowfb[k]  = 1000.*graMCNLO5fZb_noFREyLow[k];
	 graMCNLO5fZb_noFREyHighfb[k] = 1000.*graMCNLO5fZb_noFREyHigh[k];

// 	 graMCNLOYfb[k]      = 1000000.*graMCNLOY[k+offset];
// 	 graMCNLOEyLowfb[k]  = 1000000.*graMCNLOEyLow[k+offset];
// 	 graMCNLOEyHighfb[k] = 1000000.*graMCNLOEyHigh[k+offset];
	 graMCNLOYfb[k]      = 1000.*graMCNLOY[k+offset];
	 graMCNLOEyLowfb[k]  = 1000.*graMCNLOEyLow[k+offset];
	 graMCNLOEyHighfb[k] = 1000.*graMCNLOEyHigh[k+offset];
	 graMCNLOMPIYfb[k]      = 1000.*graMCNLOMPIY[k+offset];
	 graMCNLOMPIEyLowfb[k]  = 1000.*graMCNLOMPIEyLow[k+offset];
	 graMCNLOMPIEyHighfb[k] = 1000.*graMCNLOMPIEyHigh[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);
     }
     // to compute NP total syst errors 
     TGraphAsymmErrors *grMCFMNPerr  = new TGraphAsymmErrors(Npoints, 
							     grX,                             MCFMNP_Y, 
							     grEX,                            grEX,
							     MCFM_errNP_L,                    MCFM_errNP_H);
     TGraphAsymmErrors *grMCFM8NPerr  = new TGraphAsymmErrors(Npoints, 
							     grX,                             MCFM8NP_Y, 
							     grEX,                            grEX,
							     MCFM8_errNP_L,                   MCFM8_errNP_H);
     TGraphAsymmErrors *grMCFM23NPerr  = new TGraphAsymmErrors(Npoints, 
							     grX,                             MCFM23NP_Y, 
							     grEX,                            grEX,
							     MCFM23_errNP_L,                  MCFM23_errNP_H);
     // to compute total stat. err. on MCFM w/o NP 
     TGraphAsymmErrors *grMCFMstaN  = new TGraphAsymmErrors(Npoints, 
							    grX,                             grMCFM->GetY(), 
							    grEX,                            grEX,
							    grMCFMsta->GetEYlow(),           grMCFMsta->GetEYhigh());
     TGraphAsymmErrors *grMCFM8staN = new TGraphAsymmErrors(Npoints, 
							    grX,                             grMCFM8->GetY(), 
							    grEX,                            grEX,
//							    MCFM8_EYlstat,                   MCFM8_EYhstat);
							    grMCFM8sta->GetEYlow(),          grMCFM8sta->GetEYhigh());
     TGraphAsymmErrors *grMCFM23staN = new TGraphAsymmErrors(Npoints, 
							     grX,                             grMCFM23->GetY(), 
							     grEX,                            grEX,
//							     MCFM23_EYlstat,                  MCFM23_EYhstat);
							     grMCFM23sta->GetEYlow(),         grMCFM23sta->GetEYhigh());
     // MCFM with NP and stat err.s 
     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);
     // MCFM with NP and stat+alpha_s err. 
     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);
     //  MCFM with NP and stat+alpha_s+pdf err. 
     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);
     // MCFM with NP and stat+alpha_s+pdf+scale err. 
     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);
     // MCFM with NP and stat+alpha_s+pdf+scale+NPc err. 
     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* grMCFMNP_totquadrsum;
     grMCFMNP_totquadrsum       = new TGraphAsymmErrors(Npoints, 
							grX,                             MCFMNP_Y, 
							grEX,                            grEX,
							MCFMNP_EYl_totquadrsum,          MCFMNP_EYh_totquadrsum);
     TGraphAsymmErrors* grMCFM8NP_totquadrsum;
     grMCFM8NP_totquadrsum      = new TGraphAsymmErrors(Npoints, 
							grX,                             MCFM8NP_Y, 
							grEX,                            grEX,
							MCFM8NP_EYl_totquadrsum,         MCFM8NP_EYh_totquadrsum);
     TGraphAsymmErrors* grMCFM23NP_totquadrsum;
     grMCFM23NP_totquadrsum     = new TGraphAsymmErrors(Npoints, 
							grX,                             MCFM23NP_Y, 
							grEX,                            grEX,
							MCFM23NP_EYl_totquadrsum,        MCFM23NP_EYh_totquadrsum);

     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 *graMCNLO5fZj_MIfb   = new TGraphAsymmErrors(grSherpa->GetN(), 
								 grX,                    graMCNLO5fZj_MIYfb, 
								 grEX,                   grEX,
								 graMCNLO5fZj_MIEyLowfb, graMCNLO5fZj_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 ==================="<<hname<<std::endl;
     Double_t xsMCFMsta                                            = ComputeTotalCrossSection(grMCFMstaN);
     Double_t ppp   = getStatErr(grMCFMstaN);
     std::cout<<"====== MCFM MSTW2008 sta err NO NP ==================="<<hname<<std::endl;
     Double_t xsMCFM8sta                                            = ComputeTotalCrossSection(grMCFM8staN);
     std::cout<<"====== MCFM NNPDF23  sta err NO NP ==================="<<hname<<std::endl;
     Double_t xsMCFM23sta                                           = ComputeTotalCrossSection(grMCFM23staN);
     std::cout<<"======================================================"<<std::endl;
     std::cout<<"====== SHERPA no NP =================================="<<hname<<std::endl;
     Double_t xsSHERPA                                             = ComputeTotalCrossSection(grSherpafb);
     std::cout<<"======================================================"<<std::endl;
     std::cout<<"====== PYTHIA no NP =================================="<<hname<<std::endl;
     Double_t xsPYTHIA                                             = ComputeTotalCrossSection(grPythiafb);
     std::cout<<"======================================================"<<std::endl;
     std::cout<<"====== MCFM sta err =================================="<<hname<<std::endl;
     Double_t xsMCFMNPsta                                          = ComputeTotalCrossSection(grMCFMNPsta);
     //
     Double_t errStat   = getStatErr(grMCFMNPsta);
     Double_t errStat8  = getStatErr(grMCFM8NPsta);
     Double_t errStat23 = getStatErr(grMCFM23NPsta);
     for (int k=0; k<grMCFMas->GetN(); ++k) {
	 grMCFMas ->SetPointEXlow(k, grMCFMNPsta->GetEXlow()[k]);
	 grMCFMsca->SetPointEXlow(k, grMCFMNPsta->GetEXlow()[k]);
	 grMCFMpdf->SetPointEXlow(k, grMCFMNPsta->GetEXlow()[k]);
	 grMCFMas ->SetPointEXhigh(k,grMCFMNPsta->GetEXlow()[k]);
	 grMCFMsca->SetPointEXhigh(k,grMCFMNPsta->GetEXlow()[k]);
	 grMCFMpdf->SetPointEXhigh(k,grMCFMNPsta->GetEXlow()[k]);

	 grMCFM8as ->SetPointEXlow(k, grMCFMNPsta->GetEXlow()[k]);
	 grMCFM8sca->SetPointEXlow(k, grMCFMNPsta->GetEXlow()[k]);
	 grMCFM8pdf->SetPointEXlow(k, grMCFMNPsta->GetEXlow()[k]);
	 grMCFM8as ->SetPointEXhigh(k,grMCFMNPsta->GetEXlow()[k]);
	 grMCFM8sca->SetPointEXhigh(k,grMCFMNPsta->GetEXlow()[k]);
	 grMCFM8pdf->SetPointEXhigh(k,grMCFMNPsta->GetEXlow()[k]);

	 grMCFM23as ->SetPointEXlow(k, grMCFMNPsta->GetEXlow()[k]);
	 grMCFM23sca->SetPointEXlow(k, grMCFMNPsta->GetEXlow()[k]);
	 grMCFM23pdf->SetPointEXlow(k, grMCFMNPsta->GetEXlow()[k]);
	 grMCFM23as ->SetPointEXhigh(k,grMCFMNPsta->GetEXlow()[k]);
	 grMCFM23sca->SetPointEXhigh(k,grMCFMNPsta->GetEXlow()[k]);
	 grMCFM23pdf->SetPointEXhigh(k,grMCFMNPsta->GetEXlow()[k]);
     }
     Double_t errAlphasUP = getSystErr(grMCFMas, "UP");
     Double_t errAlphasDW = getSystErr(grMCFMas, "DW");
     Double_t errAlphas8UP = getSystErr(grMCFM8as, "UP");
     Double_t errAlphas8DW = getSystErr(grMCFM8as, "DW");
     Double_t errAlphas23UP = getSystErr(grMCFM23as, "UP");
     Double_t errAlphas23DW = getSystErr(grMCFM23as, "DW");
     Double_t errScaUP = getSystErr(grMCFMsca, "UP");
     Double_t errScaDW = getSystErr(grMCFMsca, "DW");
     Double_t errSca8UP = getSystErr(grMCFM8sca, "UP");
     Double_t errSca8DW = getSystErr(grMCFM8sca, "DW");
     Double_t errSca23UP = getSystErr(grMCFM23sca, "UP");
     Double_t errSca23DW = getSystErr(grMCFM23sca, "DW");
     Double_t errPdfUP = getSystErr(grMCFMpdf, "UP");
     Double_t errPdfDW = getSystErr(grMCFMpdf, "DW");
     Double_t errPdf8UP = getSystErr(grMCFM8pdf, "UP");
     Double_t errPdf8DW = getSystErr(grMCFM8pdf, "DW");
     Double_t errPdf23UP = getSystErr(grMCFM23pdf, "UP");
     Double_t errPdf23DW = getSystErr(grMCFM23pdf, "DW");
     Double_t errNPcUP = getSystErr(grMCFMNPerr, "UP");
     Double_t errNPcDW = getSystErr(grMCFMNPerr, "DW");
     Double_t errNPc8UP = getSystErr(grMCFM8NPerr, "UP");
     Double_t errNPc8DW = getSystErr(grMCFM8NPerr, "DW");
     Double_t errNPc23UP = getSystErr(grMCFM23NPerr, "UP");
     Double_t errNPc23DW = getSystErr(grMCFM23NPerr, "DW");
     //
     std::cout<<"====== MCFM sta_as err ==============================="<<hname<<std::endl;
     Double_t xsMCFMNPsta_as                                       = ComputeTotalCrossSectionWithOpt(grMCFMNPsta_as, 1);
     std::cout<<"====== MCFM sta_as_pdf err ==========================="<<hname<<std::endl;
     Double_t xsMCFMNPsta_as_pdf                                   = ComputeTotalCrossSectionWithOpt(grMCFMNPsta_as_pdf, 1);
     std::cout<<"====== MCFM sta_as_pdf_sca err ======================="<<hname<<std::endl;
     Double_t xsMCFMNPsta_as_pdf_sca                               = ComputeTotalCrossSectionWithOpt(grMCFMNPsta_as_pdf_sca, 1);
     std::cout<<"====== MCFM CT10     sta_as_pdf_sca_np err ==========="<<hname<<std::endl;
     Double_t xsMCFMNPsta_as_pdf_sca_np                          = ComputeTotalCrossSectionWithOpt(grMCFMNPsta_as_pdf_sca_np, 1);
     std::cout<<"====== MCFM CT10 sta_as_pdf_sca_np qsumerr ==========="<<hname<<std::endl;
     Double_t xsMCFMNP_qsumerr                                   = ComputeTotalCrossSectionWithOpt(grMCFMNP_totquadrsum, 1);
     std::cout<<"====== MCFM CT10 with all errors ====================="<<hname<<std::endl;
     std::cout<<"==="<<xsMCFMNPsta_as_pdf_sca_np<<" +/- "<< errStat<<"(stat) +"
	      <<errAlphasUP<<"-"<<errAlphasDW<<"(alpha_s) +"
	      <<errScaUP<<"-"<<errScaDW<<"(scales) +"
	      <<errPdfUP<<"-"<<errPdfDW<<"(pdf) +"
	      <<errNPcUP<<"-"<<errNPcDW<<"(NP_corrs.)"<<std::endl;
     double totsystDW = sqrt(errAlphasDW*errAlphasDW+errScaDW*errScaDW+errPdfDW*errPdfDW+errNPcDW+errNPcDW);
     double totsystUP = sqrt(errAlphasUP*errAlphasUP+errScaUP*errScaUP+errPdfUP*errPdfUP+errNPcUP+errNPcUP);
     std::cout<<"==="<<xsMCFMNPsta_as_pdf_sca_np<<" +/- "<< errStat<<"(stat) +"<<totsystUP<<" -"<<totsystDW<<"(syst)"<<std::endl;
     std::cout<<"====== MCFM MSTW2008 sta_as_pdf_sca_np err ==========="<<hname<<std::endl;
     Double_t xsMCFM8NPsta_as_pdf_sca_np                         = ComputeTotalCrossSectionWithOpt(grMCFM8NPsta_as_pdf_sca_np, 1);
     std::cout<<"====== MCFM MSTW2008 sta_as_pdf_sca_np qsumerr========"<<hname<<std::endl;
     Double_t xsMCFM8NP_qsumerr                                   = ComputeTotalCrossSectionWithOpt(grMCFM8NP_totquadrsum, 1);
     std::cout<<"====== MCFM MSTW2008 with all errors ================="<<hname<<std::endl;
     std::cout<<"==="<<xsMCFM8NPsta_as_pdf_sca_np<<" +/- "<< errStat8<<"(stat) +"
	      <<errAlphas8UP<<"-"<<errAlphas8DW<<"(alpha_s) +"
	      <<errSca8UP<<"-"<<errSca8DW<<"(scales) +"
	      <<errPdf8UP<<"-"<<errPdf8DW<<"(pdf) +"
	      <<errNPc8UP<<"-"<<errNPc8DW<<"(NP_corrs.) "<<std::endl;
     double totsyst8DW = sqrt(errAlphas8DW*errAlphas8DW+errSca8DW*errSca8DW+errPdf8DW*errPdf8DW+errNPc8DW+errNPc8DW);
     double totsyst8UP = sqrt(errAlphas8UP*errAlphas8UP+errSca8UP*errSca8UP+errPdf8UP*errPdf8UP+errNPc8UP+errNPc8UP);
     std::cout<<"==="<<xsMCFM8NPsta_as_pdf_sca_np<<" +/- "<< errStat8<<"(stat) +"<<totsyst8UP<<" -"<<totsyst8DW<<"(syst)"<<std::endl;
     std::cout<<"====== MCFM NNPDF23  sta_as_pdf_sca_np err ==========="<<hname<<std::endl;
     Double_t xsMCFM23NPsta_as_pdf_sca_np                        = ComputeTotalCrossSectionWithOpt(grMCFM23NPsta_as_pdf_sca_np, 1);
     std::cout<<"====== MCFM NNPDF23 sta_as_pdf_sca_np qsumerr========="<<hname<<std::endl;
     Double_t xsMCFM23NP_qsumerr                                 = ComputeTotalCrossSectionWithOpt(grMCFM23NP_totquadrsum, 1);
     std::cout<<"====== MCFM NNPDF23 with all errors =================="<<hname<<std::endl;
     std::cout<<"==="<<xsMCFM23NPsta_as_pdf_sca_np<<" +/- "<< errStat23<<"(stat) +"
	      <<errAlphas23UP<<"-"<<errAlphas23DW<<"(alpha_s) +"
	      <<errSca23UP<<"-"<<errSca23DW<<"(scales) +"
	      <<errPdf23UP<<"-"<<errPdf23DW<<"(pdf) +"
	      <<errNPc23UP<<"-"<<errNPc23DW<<"(NP_corrs.) "<<std::endl;
     double totsyst23DW = sqrt(errAlphas23DW*errAlphas23DW+errSca23DW*errSca23DW+errPdf23DW*errPdf23DW+errNPc23DW+errNPc23DW);
     double totsyst23UP = sqrt(errAlphas23UP*errAlphas23UP+errSca23UP*errSca23UP+errPdf23UP*errPdf23UP+errNPc23UP+errNPc23UP);
     std::cout<<"==="<<xsMCFM23NPsta_as_pdf_sca_np<<" +/- "<< errStat23<<"(stat) +"<<totsyst23UP<<" -"<<totsyst23DW<<"(syst)"<<std::endl;
     std::cout<<"======================================================"<<std::endl; 
     std::cout<<"====== aMC@NLO no MPI 4f UCL - v3 ===================="<<hname<<std::endl;
     Double_t xsaMCNLO4f                                         = ComputeTotalCrossSection(graMCNLO4ffb);
     std::cout<<"======================================================"<<hname<<std::endl; 
     std::cout<<"====== aMC@NLO no MPI 4f Lecce - low stat ============"<<hname<<std::endl;
     Double_t xsaMCNLO4fls                                       = ComputeTotalCrossSection(graMCNLO4flsfb);
     std::cout<<"======================================================"<<std::endl; 
     std::cout<<"====== aMC@NLO no MPI 5f Zb Lecce ===================="<<hname<<std::endl;
     Double_t xsaMCNLO5fZb                                       = ComputeTotalCrossSection(graMCNLO5fZbfb);
     std::cout<<"======================================================"<<std::endl; 
     std::cout<<"====== aMC@NLO no MPI 5f Zj Lecce ===================="<<hname<<std::endl;
     Double_t xsaMCNLO5fZj                                       = ComputeTotalCrossSection(graMCNLO5fZjfb);
     std::cout<<"======================================================"<<std::endl; 
     std::cout<<"====== aMC@NLO 5f Zb MI Lecce ========================"<<hname<<std::endl;
     Double_t xsaMCNLO5fZb_MI                                    = ComputeTotalCrossSection(graMCNLO5fZb_MIfb);
     std::cout<<"===*************** diff with MCFM_sameScaleAndPDF ===="<<std::endl;
     Double_t errStat_MCNLO  = getStatErr(graMCNLO5fZb_MIfb);	
     Double_t residual = (xsaMCNLO5fZb_MI-xsMCFM8NPsta_as_pdf_sca_np)/sqrt( errStat8*errStat8 + errStat_MCNLO*errStat_MCNLO + errNPc8UP*errNPc8UP);
     std::cout<<"===****** Resid to MCFM "<<residual<<" sigma(quadr. sum of stat+stat+NPerr)   "
		<<100.*(xsaMCNLO5fZb_MI-xsMCFM8NPsta_as_pdf_sca_np)/xsMCFM8NPsta_as_pdf_sca_np<<" %"<<std::endl;
     std::cout<<"======================================================"<<std::endl; 
     std::cout<<"====== aMC@NLO 5f Zj MI Lecce ========================"<<hname<<std::endl;
     Double_t xsaMCNLO5fZj_MI                                    = ComputeTotalCrossSection(graMCNLO5fZj_MIfb);
     std::cout<<"===*************** diff with MCFM_sameScaleAndPDF ===="<<std::endl;
     errStat_MCNLO  = getStatErr(graMCNLO5fZj_MIfb);
     residual = (xsaMCNLO5fZj_MI-xsMCFM8NPsta_as_pdf_sca_np)/sqrt( errStat8*errStat8 + errStat_MCNLO*errStat_MCNLO + errNPc8UP*errNPc8UP);
     std::cout<<"===****** Resid to MCFM "<<residual<<" sigma(quadr. sum of stat+stat+NPerr)   "
		<<100.*(xsaMCNLO5fZj_MI-xsMCFM8NPsta_as_pdf_sca_np)/xsMCFM8NPsta_as_pdf_sca_np<<" %"<<std::endl;
     std::cout<<"======================================================"<<std::endl; 
     std::cout<<"====== aMC@NLO no MPI 5f Zb noFR Lecce ==============="<<hname<<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 ============================"<<hname<<std::endl;
       xsaMCNLO                                                     = ComputeTotalCrossSection(graMCNLOfb);
       std::cout<<"====== aMC@NLO NP === UCL ============================"<<hname<<std::endl;
       xsaMCNLOMPI                                                   = ComputeTotalCrossSection(graMCNLOMPIfb);
       std::cout<<"===*************** diff with MCFM_sameScaleAndPDF ===="<<std::endl;
       errStat_MCNLO  = getStatErr(graMCNLOMPIfb);
       residual = (xsaMCNLOMPI-xsMCFM8NPsta_as_pdf_sca_np)/sqrt( errStat8*errStat8 + errStat_MCNLO*errStat_MCNLO + errNPc8UP*errNPc8UP);
       std::cout<<"===****** Resid to MCFM "<<residual<<" sigma(quadr. sum of stat+stat+NPerr)   "
		<<100.*(xsaMCNLOMPI-xsMCFM8NPsta_as_pdf_sca_np)/xsMCFM8NPsta_as_pdf_sca_np<<" %"<<std::endl;
     }
     std::cout<<"======================================================"<<std::endl;
     std::cout<<"====== ALPGEN NP ====================================="<<hname<<std::endl;
     Double_t xsALPGENNP=0.;
     if (grAlpgenNP!=NULL) xsALPGENNP = ComputeTotalCrossSection(grAlpgenNPfb);
     else std::cout<<"====== "<<std::endl;
     std::cout<<"====== ALPGEN NP ====no kfactor======================="<<hname<<std::endl;
     xsALPGENNP=0.;
     if (grAlpgenNP!=NULL) xsALPGENNP = ComputeTotalCrossSection(grAlpgenNPfb_nokf);
     else std::cout<<"====== "<<std::endl;
     std::cout<<"======================================================"<<std::endl;
     std::cout<<"====== SHERPA NP ====================================="<<hname<<std::endl;
     Double_t xsSHERPANP                                           = ComputeTotalCrossSection(grSherpaNPfb);
     std::cout<<"======================================================"<<std::endl;
     std::cout<<"====== PYTHIA NP ====================================="<<hname<<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");

	 graMCNLO5fZj_MIfb->SetFillStyle(3001);
	 graMCNLO5fZj_MIfb->SetFillColor(65);
	 graMCNLO5fZj_MIfb->SetLineColor(65);
	 graMCNLO5fZj_MIfb->SetLineStyle(1);
	 graMCNLO5fZj_MIfb->SetMarkerColor(65);
	 graMCNLO5fZj_MIfb->SetMarkerStyle(21);
	 graMCNLO5fZj_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(2);
	 graMCNLO5fZjfb->SetMarkerColor(3);
	 graMCNLO5fZjfb->SetMarkerStyle(26);
	 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(graMCNLO5fZj_MIfb,"aMC@NLO 5f Zj 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 = "v3MCFMscalesND/logY"+fname;
	     else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	     else fname = "v3MCFMscalesND/"+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 = "v3MCFMscalesND/logY"+fname;
	     else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	     else fname = "v3MCFMscalesND/"+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 = "v3MCFMscalesND/logY"+fname;
	     else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	     else fname = "v3MCFMscalesND/"+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");
	 graMCNLO5fZj_MIfb->Draw("PCSAME");

	 legend->Draw();
	 if (saveAsEps)
	   {
	     fname = "fig"+hname+"_MCNLO_4fvs5f_wMPI.eps";
	     if (setLogY) fname = "v3MCFMscalesND/logY"+fname;
	     else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	     else fname = "v3MCFMscalesND/"+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 = "v3MCFMscalesND/logY"+fname;
	     else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	     else fname = "v3MCFMscalesND/"+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);
     grMCFMstaN->SetLineColor(collMCFM);
     grMCFMstaN->SetMarkerColor(colmMCFM);
     grMCFMstaN->SetLineWidth(1);
     grMCFMstaN->GetXaxis()->SetTitle(xTitle);
     grMCFMstaN->GetYaxis()->SetTitle(yTitle);

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

     grMCFM23staN->SetLineColor(collMCFM23);
     grMCFM23staN->SetMarkerColor(colmMCFM23);
     grMCFM23staN->SetLineWidth(1);
     grMCFM23staN->GetXaxis()->SetTitle(xTitle);
     grMCFM23staN->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") grMCFMstaN->GetYaxis()->SetRangeUser(0.01, 500.);

     if (hname=="j_2bjet_Mbb_binned") {
       grSherpafb->Draw("PCSAME");
       grMCFMstaN->Draw("PCSAME");
     }
     else {
       grMCFMstaN->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(grMCFMstaN, "MCFM","PE");	
     legend->Draw();

     if (saveAsEps)
       {
	 std::string fname = "fig"+hname+"_MCFM_Sherpa_Pythia_noNP.eps";
	 if (setLogY) fname = "v3MCFMscalesND/logY"+fname;
	 else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	 else fname = "v3MCFMscalesND/"+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(2);
	 grAlpgenNPfb->SetLineColor(colAlpgen);
	 grAlpgenNPfb->SetMarkerColor(colAlpgen);
	 grAlpgenNPfb->Draw("PCSAME");
	 grAlpgenNPfb_nokf->SetLineStyle(1);
	 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 = "v3MCFMscalesND/logY"+fname;
	 else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	 else fname = "v3MCFMscalesND/fig"+hname+"_aMCNLO_Sherpa_Pythia.eps";
	 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");
     }     
     graMCNLOMPIfb->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");
     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 = "v3MCFMscalesND/logY"+fname;
	 else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	 else fname = "v3MCFMscalesND/fig"+hname+"_aMCNLO45f_Sherpa_Alpgen.eps";
	 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");
     grMCFM8NP_totquadrsum->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(grMCFM8NP_totquadrsum, "MCFM mstw8 - total err", "PEL");

     legend->Draw();
     
     if (saveAsEps)
       {
	 fname = "fig"+hname+"_aMCNLO45f_Sherpa_Py_Alpgen.eps";
	 if (setLogY) fname = "v3MCFMscalesND/logY"+fname;
	 else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	 else fname = "v3MCFMscalesND/fig"+hname+"_aMCNLO45f_Sherpa_Py_Alpgen.eps";
	 cPlotSPA->SaveAs(fname.c_str());
       }


     //if (cPlotNP) delete cPlotNP;
     TCanvas *cPlotMCFMCT10NP = new TCanvas("MCFMwNPCT10","MCFMwNP_CT10");
     if (setLogY) cPlotMCFMCT10NP->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");
     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->Draw();
     
     if (saveAsEps)
       {
	 fname = "fig"+hname+"_MCFM_CT10.eps";
	 if (setLogY) fname = "v3MCFMscalesND/logY"+fname;
	 else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	 else fname = "v3MCFMscalesND/"+fname;
	 cPlotMCFMCT10NP->SaveAs(fname.c_str());
       }

     TCanvas *cPlotNP = new TCanvas("MCFM_AMC@NLO_PYTHIA_SHERPA","MCFM_AMC@NLO_PYTHIA_SHERPA");
     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 = "v3MCFMscalesND/logY"+fname;
	 else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	 else fname = "v3MCFMscalesND/fig"+hname+"_MCFM_aMCNLO_Sherpa_Pythia.eps";
	 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 = "v3MCFMscalesND/logY"+fname;
	 else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	 else fname = "v3MCFMscalesND/fig"+hname+"_MCFM_aMCNLO_Sherpa_Pythia_WITHRATIO.eps";
	 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");
//      // NNPDF23 central value with NP 
//      grMCFM23NPsta->Draw("PESAME");


     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();
     grMCFMstaN-> SetLineStyle(1);
     

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

	      grMCFMNP_totquadrsum->Write("grMCFMNP_totquadrsum");
	      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");
	      grMCFM8NP_totquadrsum->Write("grMCFM8NP_totquadrsum");
	      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");
	      grMCFM23NP_totquadrsum->Write("grMCFM23NP_totquadrsum");
	      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 = "v3MCFMscalesND/logY"+fname;
	 else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	 else fname = "v3MCFMscalesND/fig"+hname+"_MCFM2pdf.eps";
	 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");
//      // MSTW2008 central value with NP 
//      grMCFM8NPsta->Draw("PESAME");


     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();
     grMCFMstaN-> SetLineStyle(1);
     

     if (saveAsEps)
       {
	 fname = "fig"+hname+"_MCFM23pdf.eps";
	 if (setLogY) fname = "v3MCFMscalesND/logY"+fname;
	 else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	 else fname = "v3MCFMscalesND/fig"+hname+"_MCFM23pdf.eps";
	 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");
     grMCFMNP_totquadrsum->SetFillColor(5);
     //grMCFMNP_totquadrsum->Draw("E2SAME");
     grMCFMNP_totquadrsum->Draw("E2");
     //grMCFMsta->Draw("CSAME");
     grMCFM8NP_totquadrsum->SetFillStyle(3004);
     grMCFM8NP_totquadrsum->SetFillColor(collMCFM8);
     grMCFM8NP_totquadrsum->Draw("E2SAME");
     //grMCFM8sta->Draw("CSAME");
     grMCFM23NP_totquadrsum->SetFillStyle(3005);
     grMCFM23NP_totquadrsum->SetFillColor(collMCFM23);
     grMCFM23NP_totquadrsum->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(grMCFMNP_totquadrsum,   "MCFM (ct10) x NP",    "F");
     legend->AddEntry(grMCFM8NP_totquadrsum,  "MCFM (mstw2008) x NP","F");
     legend->AddEntry(grMCFM23NP_totquadrsum, "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 = "v3MCFMscalesND/logY"+fname;
	 else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	 else fname = "v3MCFMscalesND/fig"+hname+"_MCFM_w_wo_NP.eps";
	 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");
	 grMCFMNP_totquadrsum->SetFillColor(5);
	 //grMCFMNP_totquadrsum->Draw("E2SAME");
	 grMCFMNP_totquadrsum->Draw("E2");
	 //grMCFMsta->Draw("CSAME");
	 grMCFM8NP_totquadrsum->SetFillStyle(3004);
	 grMCFM8NP_totquadrsum->SetFillColor(collMCFM8);
	 grMCFM8NP_totquadrsum->Draw("E2SAME");
	 //grMCFM8sta->Draw("CSAME");
	 grMCFM23NP_totquadrsum->SetFillStyle(3005);
	 grMCFM23NP_totquadrsum->SetFillColor(collMCFM23);
	 grMCFM23NP_totquadrsum->Draw("E2SAME");
	 //grMCFM23sta->Draw("CSAME");

	 if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))   
	   {     
	     std::cout<<" DO I ENTER HERE "<<std::endl;
	     /// test ---------------{
	     graMCNLOMPIfb->SetLineColor(2);
	     graMCNLOMPIfb->SetMarkerColor(2);
	     graMCNLOMPIfb->Draw("PSAME");
	     //	     graMCNLOMPIfb->Draw("PESAME");
	     /// 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);

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

	 grdataMBayesSyst->SetMarkerColor(1);
	 grdataMBayesSyst->SetMarkerStyle(20);
	 grdataMBayesSyst->SetMarkerSize(0.1);
	 grdataMBayesTot->SetMarkerColor(1);
	 grdataMBayesTot->SetMarkerStyle(20);
	 grdataMBayesTot->SetMarkerSize(0.1);
    
	 //	 grdataMBayesTot->SetLineWidth(1.);
	 //	 grdataMBayesSyst->SetLineWidth(1.);

	 grdataMbbb->Draw("PESAME");
	 grdataMBayes->Draw("PESAME");     
	 grdataNbbb->Draw("PESAME");
	 grdataNBayes->Draw("PESAME");
	 legend1->Clear();
	 legend1->AddEntry(grMCFMNP_totquadrsum,"MCFM (ct10) x NP","F");
	 //     legend->AddEntry(grMCFMsta           ,     "MCFM (ct10) stat. err.","L");
	 legend1->AddEntry(grMCFM8NP_totquadrsum,"MCFM (mstw2008) x NP","F");
	 //     legend->AddEntry(grMCFM8sta           ,     "MCFM (mstw2008) stat. err.","L");
	 legend1->AddEntry(grMCFM23NP_totquadrsum,"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 = "v3MCFMscalesND/logY"+fname;
	     else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	     else         fname = "v3MCFMscalesND/"+fname;
	     cPlotNP_01data->SaveAs(fname.c_str());
	   }
	 //	 return;


	 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();
	 if (setLogX) c1DandR->SetLogx();
	 c1DandR->cd();
	 hemptyCopy->Draw();
	 
	 if (!skipMCFM)
	 {
	     
	     // 	 grMCFMNP_totquadrsum->Draw("E2");
	     // 	 grMCFM8NP_totquadrsum->Draw("E2SAME");
	     // 	 grMCFM23NP_totquadrsum->Draw("E2SAME");
	     grMCFM8NP_totquadrsum ->SetLineColor(4);	 
	     grMCFMNP_totquadrsum  ->SetLineColor(3);	 
	     grMCFM23NP_totquadrsum->SetLineColor(7);	 

	     grMCFM8NP_totquadrsum ->SetLineWidth(2);	 
	     grMCFMNP_totquadrsum  ->SetLineWidth(2);	 
	     grMCFM23NP_totquadrsum->SetLineWidth(2);	 
	     //
	     grMCFM8NP_totquadrsum ->SetMarkerColor(4);	 
	     grMCFMNP_totquadrsum  ->SetMarkerColor(3);	 
	     grMCFM23NP_totquadrsum->SetMarkerColor(7);	 
	     //
	     grMCFM8NP_totquadrsum ->SetMarkerSize(0.01);	 
	     grMCFMNP_totquadrsum  ->SetMarkerSize(0.01);	 
	     grMCFM23NP_totquadrsum->SetMarkerSize(0.01);	 
	     //
	     grMCFM8NP_totquadrsum ->SetFillColor(4);
	     grMCFM8NP_totquadrsum ->SetFillStyle(3001);
	     //grMCFMNP_totquadrsum  ->Draw("CX");
	     //grMCFM8NP_totquadrsum ->Draw("CX");//tornaqui
	     grMCFM8NPsta->SetFillColor(4);
	     grMCFM8NPsta->SetFillStyle(3001);
	     grMCFM8NPsta->SetMarkerSize(grMCFM8NP_totquadrsum ->GetMarkerSize());	 
	     grMCFM8NPsta->SetMarkerColor(grMCFM8NP_totquadrsum ->GetMarkerColor());	 
	     grMCFM8NPsta->SetLineWidth(grMCFM8NP_totquadrsum ->GetLineWidth());	 
	     grMCFM8NPsta->SetLineColor(grMCFM8NP_totquadrsum ->GetLineColor());	 
	     grMCFM8NPsta->Draw("E2");
	     grMCFM8NPsta->Draw("EPZ");
	     if (showLines) grMCFM8NPsta->Draw("CX");
       }
     
	 grSherpaNPfb ->SetLineWidth(3);
	 grSherpaNPfb ->SetMarkerSize(0.01);
	 //	 grSherpaNPfb ->Draw("CX");	
	 grSherpaNPfb->SetFillColor(3);
	 grSherpaNPfb->SetFillStyle(3001);
	 grSherpaNPfb ->Draw("E2");
	 grSherpaNPfb ->Draw("EPZ");
	 if (showLines) grSherpaNPfb ->Draw("CX");
	 //
	 if (grAlpgenNP!=NULL)
	 {
	     grAlpgenNPfb_nokf ->SetLineWidth(3);
	     grAlpgenNPfb_nokf ->SetMarkerSize(0.01);
	     grAlpgenNPfb_nokf->SetFillColor(grAlpgenNPfb_nokf->GetLineColor());
	     grAlpgenNPfb_nokf->SetFillStyle(3001);
	 //	     grAlpgenNPfb_nokf ->Draw("CX");
	     grAlpgenNPfb_nokf ->Draw("E2");
	     grAlpgenNPfb_nokf ->Draw("EPZ");
	     if (showLines) grAlpgenNPfb_nokf ->Draw("CX");
	 }
	 
	 //grMCFM8NP_totquadrsum ->Draw("E4");
	 //grMCFM23NP_totquadrsum->Draw("CX");

	 if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))   
	   {     
	     /// test ---------------{
	     graMCNLOMPIfb->SetMarkerSize(0.01);
	     graMCNLOMPIfb->SetLineWidth(2);
	     graMCNLOMPIfb->SetLineStyle(1);
	     graMCNLOMPIfb->SetLineColor(2);
	     graMCNLOMPIfb->SetFillColor(2);
	     graMCNLOMPIfb->SetFillStyle(3001);
	     graMCNLOMPIfb->SetMarkerColor(2);
	     //	     graMCNLOMPIfb->TGraphPainter::Draw("CX");
	     graMCNLOMPIfb->TGraphPainter::Draw("E2");
	     graMCNLOMPIfb->TGraphPainter::Draw("EPZ");
	     if (showLines) graMCNLOMPIfb->TGraphPainter::Draw("CX");
	   }
	 graMCNLO5fZb_MIfb->SetLineWidth(2.);
	 graMCNLO5fZb_MIfb->SetMarkerSize(0.1);
	 if (showZbInPlots) graMCNLO5fZb_MIfb->Draw("PCSAME");
	 graMCNLO5fZj_MIfb->SetLineWidth(2.);
	 graMCNLO5fZj_MIfb->SetMarkerSize(0.3);
	 graMCNLO5fZj_MIfb->Draw("PE2");
	 graMCNLO5fZj_MIfb->Draw("EPZ");
	 if (showLines) graMCNLO5fZj_MIfb->Draw("CX");

	 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(26);
	 grdataNbbb->SetMarkerSize(1.5);
	 grdataNBayes->SetMarkerStyle(23);
	 grdataMbbb->SetMarkerStyle(24);
	 grdataMbbb->SetMarkerSize(1.5);
	 grdataMBayesTot->SetMarkerSize(1.);
	 grdataMBayes->SetMarkerStyle(20);
    	
	 //gStyle->SetErrorX(0);// controls the size of the error along x. dx = 0 removes the error along x.
	 // gStyle->SetEndErrorSize(10.);

	 //	 if (Mbbb)   grdataMbbb->Draw("PESAME");
	 if (MBayes) {
// 	   grdataMBayesTot->SetLineWidth(2);
// 	   grdataMBayesSyst->SetLineWidth(2);
//	   grdataMBayesTot->Draw("PESAMEZ");     
//	   grdataMBayesSyst->Draw("PESAME");   
	 }
	 if (Nbbb)   grdataNbbb->Draw("PESAME");
	 if (NBayes) {
	     grdataNBayesTot->Draw("PESAME");
	     //grdataNBayesSyst->Draw("PESAME");
	 }
	 //gStyle->SetErrorX(1);// controls the size of the error along x. dx = 0 removes the error along x.
	 //gStyle->SetEndErrorSize(0.);


	 
	 legendDMC->Clear();
 	 if (!skipMCFM) legendDMC->AddEntry(grMCFM8NP_totquadrsum,"MCFM (mstw2008) x NP","F");
	 //legendDMC->AddEntry(grMCFMNP_totquadrsum,"MCFM (ct10) x NP","L");
	 //legendDMC->AddEntry(grMCFM23NP_totquadrsum,"MCFM (nnpdf23) x NP","L");
	 if (hnameaMCNLO != "Z_y" || (!oldaMCNLO))  
	   {
	     //legendDMC->AddEntry(graMCNLOMPIfb,            "aMC@NLO Zbb, 4fns","L");
	     legendDMC->AddEntry(graMCNLOMPIfb,            "aMC@NLO Zbb, 4fns","F");
	     if (showMCNLOnoMPI) legendDMC->AddEntry(graMCNLOfb,               "aMC@NLO no MPI","L");
	   }
	 if (showZbInPlots) legendDMC->AddEntry(graMCNLO5fZb_MIfb,"aMC@NLO Zb, 5fns","L");
	 //	 legendDMC->AddEntry(graMCNLO5fZj_MIfb,"aMC@NLO Zj, 5fns","L");
	 legendDMC->AddEntry(graMCNLO5fZj_MIfb,"aMC@NLO Zj, 5fns","F");
	 legendDMC->AddEntry(grSherpaNPfb,"SHERPA","L");
	 if (grAlpgenNP!=NULL) legendDMC->AddEntry(grAlpgenNPfb_nokf,"Alpgen","L");
	 
	 if (MBayes) legendDMC->AddEntry(grdataMBayesTot,          "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");	     
	     myText(0.65,0.07,1,"Z p_{T} > 20 GeV");	     
	 }
	 else if (hname=="j_bjet_ZdR620_binned") 
	 {
	     ATLAS_LABEL(0.2,0.85);
	     myText(0.33,0.85,1,"Internal");
	     myText(0.65,0.85,1,"Z+b selection");	     
	     myText(0.65,0.77,1,"Z p_{T} > 20 GeV");	     
	 }
	 else if (hname=="j_bjet_Pt_binned")
	 {
	     ATLAS_LABEL(0.2,0.87);
	     myText(0.33,0.87,1,"Internal");
	     myText(0.6,0.5,1,"Z+b selection");
	 }
	 else 
	 {
	     ATLAS_LABEL(0.2,0.87);
	     myText(0.33,0.87,1,"Internal");
	     myText(0.2,0.07,1,"Z+b selection");
	     if (hname=="j_bjet_DY20_binned") myText(0.2,0.15,1,"Z p_{T} > 20 GeV");	     
	 }
	 
	 

	 c2DandR->cd();
	 
	 // syst only - assuming 15% 
	 TGraphAsymmErrors* grdataSyst;
	 // stat + syst		      
	 TGraphAsymmErrors* grdataStat;

	 if (MBayes) {
	     grdataSyst = grdataMBayesSyst;	     
	     grdataStat = grdataMBayesTot;	     
	 }
	 if (NBayes) {
	     grdataSyst = grdataNBayesSyst;	     
	     grdataStat = grdataNBayesTot;	     
	 }
	 
	 grdataSyst->SetFillColor(1);
	 grdataSyst->SetFillStyle(3005);
	 grdataStat->SetFillColor(5);
	 

	 if (setLogX) c2DandR->SetLogx();
	 hemptyRat->Draw();

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


	     grMCFM8NP_totquadrsum->SetLineColor(4);	 
	     grMCFMNP_totquadrsum->SetLineColor(3);	 
	     grMCFM23NP_totquadrsum->SetLineColor(7);	 
	     //
	     grMCFM8NP_totquadrsum ->SetLineWidth(2);	 
	     grMCFMNP_totquadrsum  ->SetLineWidth(4);	 
	     grMCFM23NP_totquadrsum->SetLineWidth(4);	 
	     grMCFMNP_totquadrsum  ->SetLineStyle(2);	 
	     grMCFM23NP_totquadrsum->SetLineStyle(9);	 
	     //
	     grMCFM8NP_totquadrsum->SetMarkerColor(4);	 
	     grMCFMNP_totquadrsum->SetMarkerColor(3);	 
	     grMCFM23NP_totquadrsum->SetMarkerColor(9);	 
	     //
	     grMCFM8NP_totquadrsum->SetMarkerSize(0.01);	 
	     grMCFMNP_totquadrsum->SetMarkerSize(0.01);	 
	     grMCFM23NP_totquadrsum->SetMarkerSize(0.01);	 

	 if (MBayes) {
	   rat0D = getGraphRatio(grdataMBayes,grdataStat);
	   rat1D = getGraphRatio(grdataMBayes,grdataSyst);
	   if (!skipMCFM)
	   {
//   	       grMCFMNP_totquadrsum->SetLineStyle(3);
//  	       grMCFM23NP_totquadrsum->SetLineStyle(7);
//  	       grMCFMNP_totquadrsum->SetLineWidth(4);
//  	       grMCFM23NP_totquadrsum->SetLineWidth(4);
	       
	       rat2D = getGraphRatio(grdataMBayes,grMCFMNP_totquadrsum);
	       rat3D = getGraphRatio(grdataMBayes,grMCFM8NP_totquadrsum);
	       rat4D = getGraphRatio(grdataMBayes,grMCFM23NP_totquadrsum);
	   }
	   rat5D = getGraphRatio(grdataMBayes,graMCNLOMPIfb);
	   rat6D = getGraphRatio(grdataMBayes,grSherpaNPfb);
	   rat8D = getGraphRatio(grdataMBayes,graMCNLO5fZb_MIfb);
	   rat9D = getGraphRatio(grdataMBayes,graMCNLO5fZj_MIfb);
	   if (grAlpgenNP!=NULL) rat7D = getGraphRatio(grdataMBayes,grAlpgenNPfb_nokf);
	 }

	 rat0D->Draw("E2");
	 rat1D->Draw("E2SAME");
	 //rat2D->Draw("EPSAME");
	 /// 5 and 9 are aMC@NLO 4f and 5f 
	 //rat5D->Draw("EPSAME");
	 rat5D->Draw("E2");
	 rat5D->Draw("EPZ");
	 rat6D->Draw("EPSAMEZ");
	 rat7D->Draw("EPSAMEZ");
	 if (showZbInPlots) rat8D->Draw("EPSAMEZ");
	 //rat9D->Draw("EPSAME");
 	 rat9D->Draw("E2");
 	 rat9D->Draw("EPZ");
	 if (!skipMCFM) rat3D->Draw("EPSAMEZ");
	 //rat4D->Draw("EPSAME");
	 legendRATD->Clear();
 	 legendRATD->AddEntry(rat0D,"Data stat+syst err.","F");
 	 legendRATD->AddEntry(rat1D,"Data syst err.","F");
	 legendRATS->Clear();
	 if (!skipMCFM) legendRATS->AddEntry(rat3D,"MCFM(mstw8)xNP, stat+syst err.","L");
	 legendRATS->AddEntry(rat5D,"aMC@NLO Zbb, 4fns","F");
	 if (showZbInPlots) legendRATS->AddEntry(rat8D,"aMC@NLO Zb, 5fns","L");
	 legendRATS->AddEntry(rat9D,"aMC@NLO Zj, 5fns","F");
	 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 = "v3MCFMscalesND/logY"+fname;
	     else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	     else fname = "v3MCFMscalesND/"+fname;
	     cDandR->SaveAs(fname.c_str());
	   }


	 ///////////////////////////// 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 }


	 if (!skipMCFM)
	 {
	     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();
	     if (setLogX) {
		 c1DandRPDF->SetLogx();
		 c2DandRPDF->SetLogx();
		 c3DandRPDF->SetLogx();
	     }
	     hemptyCopy->Draw();

	     //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 ->SetLineColor(0);	 
	     grMCFMNPsta_as_pdf_sca_np  ->SetLineColor(0);	 
	     grMCFM23NPsta_as_pdf_sca_np->SetLineColor(0);	 

	     //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 ->SetLineWidth(0.);	 
	     grMCFMNPsta_as_pdf_sca_np  ->SetLineWidth(0.);	 
	     grMCFM23NPsta_as_pdf_sca_np->SetLineWidth(0.);	 
	     
	     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);	 
	     //
	     //
	     grMCFMNP_totquadrsum  ->Draw("CX");
	     grMCFM23NP_totquadrsum->Draw("CX");
	     grMCFM8NP_totquadrsum ->Draw("CX");

	     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(20);
    

	     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(grMCFM8NP_totquadrsum,"MCFM (mstw2008) x NP","L");
	     legendDMC->AddEntry(grMCFMNP_totquadrsum,"MCFM (ct10) x NP","L");
	     legendDMC->AddEntry(grMCFM23NP_totquadrsum,"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(grdataMBayesTot,             "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.85);
		 myText(0.33,0.85,1,"Internal");
		 myText(0.65,0.85,1,"Z+b selection");	     
	     }
	     else if (hname=="j_bjet_Pt_binned")
	     {
		 ATLAS_LABEL(0.2,0.87);
		 myText(0.33,0.87,1,"Internal");
		 myText(0.6,0.5,1,"Z+b selection");
	     }
	     else 
	     {
		 ATLAS_LABEL(0.2,0.87);
		 myText(0.33,0.87,1,"Internal");
		 myText(0.2,0.07,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->SetFillColor(kWhite);
	     grMCFM8NPsta->SetFillStyle(1001);
	     grMCFM8NPsta->SetLineColor(4);
	     grMCFM8NPsta->SetLineWidth(2);
	     TGraphAsymmErrors* rat3Dstat = getGraphRatio(grMCFM8NPsta,grMCFM8NPsta);
	     TGraphAsymmErrors* rat3Dstat_as = getGraphRatio(grMCFM8NPsta,grMCFM8NPsta_as);
	     TGraphAsymmErrors* rat3Dstat_as_pdf = getGraphRatio(grMCFM8NPsta,grMCFM8NPsta_as_pdf);
	     TGraphAsymmErrors* rat3Dstat_as_pdf_sca = getGraphRatio(grMCFM8NPsta,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(grMCFM8NPsta,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);
	     /////////////////////// no need here /////////// 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");
	     rat0D->Draw("PE");


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

	     if (legendRATDpdf==0) {
		 //		 if( "j_bjet_Pt_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 (legendRATD1pdf==0) {
		 //	 if( "j_bjet_Pt_binned"==hname) legendRATD=new TLegend(0.175,0.55,0.475,0.90);
		 legendRATD1pdf=new TLegend(0.17,0.80,0.5,0.97);  
	     }
	     legendRATD1pdf->SetShadowColor(1);
	     legendRATD1pdf->SetBorderSize(0);
	     legendRATD1pdf->SetFillStyle(0);
	     //legend->SetFillColor(kWhite);
	     legendRATD1pdf->SetTextSize(0.06);  

	     if (legendRATSpdf==0) {
		 double dx = 0.48;
		 double dy = 0.32;
		 
		 //legendRATS=new TLegend(0.5,0.70,0.98,0.97);   // with 4 theories 
		 if ("j_bjet_Pt_binned"==hname || "j_bjet_Y_binned"==hname ) 
		     legendRATSpdf=new TLegend(0.3,0.03,0.3+dx,0.03+dy);
		 else if ( "j_bjet_DY20_binned"==hname ||  "j_bjet_Yboost_binned"==hname 
			   || "j_bjet_ZY_binned"==hname ) 
		     legendRATSpdf=new TLegend(0.2,0.03,0.2+dx,0.03+dy);
		 else if ( "j_bjet_ZPt_binned"==hname ) 
		     legendRATSpdf=new TLegend(0.15,0.03,0.15+dx,0.03+dy);
		 else if ( "j_bjet_ZdR620_binned"==hname ) 
		     legendRATSpdf=new TLegend(0.15,0.18,0.15+dx,0.18+dy);
		 else legendRATSpdf=new TLegend(0.5,0.65,0.5+dx,0.65+dy);   // with 2 theories 
	     }
	     if (legendRATSpdf1==0) {
		 double dx = 0.48;
		 double dy = 0.20;
		 //legendRATS=new TLegend(0.5,0.70,0.98,0.97);   // with 4 theories 
		 if ("j_bjet_Pt_binned"==hname || "j_bjet_Y_binned"==hname ) 
		     legendRATSpdf1=new TLegend(0.3,0.23,0.3+dx,0.23+dy);
		 else if ( "j_bjet_DY20_binned"==hname || "j_bjet_Yboost_binned"==hname 
			   || "j_bjet_ZY_binned"==hname ) 
		     legendRATSpdf1=new TLegend(0.2,0.23,0.2+dx,0.23+dy);
		 else if ( "j_bjet_ZPt_binned"==hname ) 
		     legendRATSpdf1=new TLegend(0.15,0.23,0.15+dx,0.23+dy);
		 else if ( "j_bjet_ZdR620_binned"==hname ) 
		     legendRATSpdf1=new TLegend(0.15,0.43,0.15+dx,0.43+dy);
		 else legendRATSpdf1=new TLegend(0.5,0.75,0.5+dx,0.75+dy);   // with 2 theories 
	     }
	     legendRATSpdf->SetShadowColor(1);
	     legendRATSpdf->SetBorderSize(0);
	     legendRATSpdf->SetFillStyle(0);
	     legendRATSpdf->SetTextSize(0.06);  
	     legendRATSpdf1->SetShadowColor(1);
	     legendRATSpdf1->SetBorderSize(0);
	     legendRATSpdf1->SetFillStyle(0);
	     legendRATSpdf1->SetTextSize(0.06);  



	     legendRATDpdf->Clear();
	     legendRATDpdf->AddEntry(rat0D,"Data stat+syst err.","PEL");
	     ////////////////////// not needed here ////////////    legendRATDpdf->AddEntry(rat1D,"Data syst err.","F");
	     legendRATSpdf->Clear();
	     legendRATSpdf->AddEntry(rat3Dstat,                "MCFM(mstw8)xNP, stat. err.","PEL");
	     legendRATSpdf->AddEntry(rat3Dstat_as,             "MCFM(mstw8)xNP, #alpha_{s} err.","F");
	     legendRATSpdf->AddEntry(rat3Dstat_as_pdf,         "MCFM(mstw8)xNP, pdf err.","F");
	     legendRATSpdf->AddEntry(rat3Dstat_as_pdf_sca,     "MCFM(mstw8)xNP, scale err.","F");
	     legendRATSpdf->AddEntry(rat3Dstat_as_pdf_sca_np,  "MCFM(mstw8)xNP, NP err.","F");
	     legendRATDpdf->Draw();
	     legendRATSpdf->Draw();


	     c3DandRPDF->cd();

	     legendRATD1pdf->Clear();
	     legendRATD1pdf->AddEntry(rat0D,"Data stat+syst err.","F");
	     legendRATD1pdf->AddEntry(rat1D,"Data syst err.","F");

	     hemptyRat->Draw();
	     rat0D->Draw("E2");
	     rat1D->Draw("E2SAME");

// 	       grMCFMNP_totquadrsum->SetLineStyle(3);
// 	       grMCFM23NP_totquadrsum->SetLineStyle(7);
// 	       grMCFMNP_totquadrsum->SetLineWidth(4);
// 	       grMCFM23NP_totquadrsum->SetLineWidth(4);

//	     rat2D->Draw("CSAME");
//	     rat4D->Draw("CSAME");
	     TGraphAsymmErrors* rat2Dn= rat2D->Clone();
	     TGraphAsymmErrors* rat4Dn= rat4D->Clone();
	     TGraphAsymmErrors* rat3Dn= rat3D->Clone();
	     rat2Dn->SetLineStyle(1);
	     rat4Dn->SetLineStyle(1);
	     rat2Dn->SetLineWidth(2);
	     rat4Dn->SetLineWidth(2);
	     rat2Dn->SetMarkerColor(rat2Dn->GetLineColor());
	     rat4Dn->SetMarkerColor(rat4Dn->GetLineColor());
	     rat2Dn->SetMarkerSize(1.5);
	     rat4Dn->SetMarkerSize(1.5);
	     rat2Dn->SetMarkerStyle(20);
	     rat4Dn->SetMarkerStyle(23);
	     rat2Dn->Draw("EP");
	     rat4Dn->Draw("EP");

	     rat3Dn->Draw("EPCSAME0");

	     legendRATSpdf1->Clear();
	     legendRATSpdf1->AddEntry(rat4Dn,"MCFM(nnpdf23)xNP","P");
	     legendRATSpdf1->AddEntry(rat3Dn,"MCFM(mstw8)xNP","LP");
	     legendRATSpdf1->AddEntry(rat2Dn,"MCFM(ct10)xNP","P");
	     legendRATD1pdf->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 = "v3MCFMscalesND/logY"+fname;
		 else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
		 else fname = "v3MCFMscalesND/"+fname;
		 cDandRPDF->SaveAs(fname.c_str());
	     }
	 }
	 
	 

	 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 = "v3MCFMscalesND/logY"+fname;
	     else if (setLogX) fname = "v3MCFMscalesND/logX"+fname;
	     else fname = "v3MCFMscalesND/"+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
*///////////////