C++ RooRealVar::getVal方法代码示例

//____________________________________
void DoSPlot(RooWorkspace* ws){
  std::cout << "Calculate sWeights" << std::endl;

  // get what we need out of the workspace to do the fit
  RooAbsPdf* model = ws->pdf("model");
  RooRealVar* zYield = ws->var("zYield");
  RooRealVar* qcdYield = ws->var("qcdYield");
  RooDataSet* data = (RooDataSet*) ws->data("data");

  // fit the model to the data.
  model->fitTo(*data, Extended() );

  // The sPlot technique requires that we fix the parameters
  // of the model that are not yields after doing the fit.
  RooRealVar* sigmaZ = ws->var("sigmaZ");
  RooRealVar* qcdMassDecayConst = ws->var("qcdMassDecayConst");  
  sigmaZ->setConstant();
  qcdMassDecayConst->setConstant();


  RooMsgService::instance().setSilentMode(true);


  // Now we use the SPlot class to add SWeights to our data set
  // based on our model and our yield variables
  RooStats::SPlot* sData = new RooStats::SPlot("sData","An SPlot",
		            *data, model, RooArgList(*zYield,*qcdYield) );


  // Check that our weights have the desired properties

  std::cout << "Check SWeights:" << std::endl;


  std::cout << std::endl <<  "Yield of Z is " 
	    << zYield->getVal() << ".  From sWeights it is "
	    << sData->GetYieldFromSWeight("zYield") << std::endl;


  std::cout << "Yield of QCD is " 
	    << qcdYield->getVal() << ".  From sWeights it is "
	    << sData->GetYieldFromSWeight("qcdYield") << std::endl
	    << std::endl;

  for(Int_t i=0; i < 10; i++)
    {
      std::cout << "z Weight   " << sData->GetSWeight(i,"zYield") 
		<< "   qcd Weight   " << sData->GetSWeight(i,"qcdYield") 
		<< "  Total Weight   " << sData->GetSumOfEventSWeight(i) 
		<< std::endl;
    }

  std::cout << std::endl;

  // import this new dataset with sWeights
 std::cout << "import new dataset with sWeights" << std::endl;
 ws->import(*data, Rename("dataWithSWeights"));


}

void PEs(RooAbsPdf *iGen,RooAbsPdf *iFit,int iN,int iNEvents,RooRealVar &iVar,RooRealVar &iSig,RooRealVar &iMean,
	 RooRealVar &iScale,RooRealVar &iRes) { 
  double iM0 = iMean.getVal(); double iS0 = iSig.getVal();
  //iScale.setVal(iMeanScale); iRes.setVal(iSigScale);
  TRandom1 *lRand = new TRandom1(0xDEADBEEF);
  TNtuple * lDN= new TNtuple( "xxx","xxx","ntot:m_r:m:merr:sig_r:sig:sigerr");
  for(int i0=0;i0<iN;i0++){
    if(i0 % 10 == 0) cout << "+++++++++++++++++++++++++++ running ======> " << i0 << endl;
    int lN    = lRand->Poisson(iNEvents);
    RooDataSet * lSignal  = iGen->generate(iVar,lN);
    iMean.setVal(iM0); iSig.setVal(iS0);
    iFit->fitTo(*lSignal,Strategy(1));//,Save(kTRUE),PrintLevel(1));
    if(iMean.getError() < 0.05) iFit->fitTo(*lSignal,Strategy(2));
    Float_t values[]={
      (Float_t) lN,
      (Float_t) 90.78/iScale.getVal(),
      (Float_t) iMean.getVal(),
      (Float_t) iMean.getError(),
      (Float_t) iSig.getVal(),
      (Float_t) iSig.getVal(),
      (Float_t) iSig.getError()
    };
    lDN->Fill(values);
  }
  TFile *lF = new TFile("XXX.root","RECREATE");
  lDN->Write();
  lF->Close();
}

TGraph *graphLH(std::string nuisname, double err ){

	w->loadSnapshot("bestfitall"); // SetTo BestFit values as start

	// Get The parameter we want 
	RooRealVar *nuis =(RooRealVar*) w->var(nuisname.c_str());
	double bf = nuis->getVal();
	double nll_0=nll->getVal();


	TGraph *gr = new TGraph(2*npoints+1);
	for (int i=-1*npoints;i<=npoints;i++){
		nuis->setVal(bf+err*( ((float)i)*nsigma/npoints));
		double nll_v = nll->getVal();
		gr->SetPoint(i+npoints,nuis->getVal(),nll_v-nll_0);
	}

	gr->SetTitle("");
	gr->GetYaxis()->SetTitle("NLL - obs data");
	gr->GetYaxis()->SetTitleOffset(1.1);
	gr->GetXaxis()->SetTitleSize(0.05);
	gr->GetYaxis()->SetTitleSize(0.05);
	gr->GetXaxis()->SetTitle(nuisname.c_str());
	gr->SetLineColor(4);
	gr->SetLineWidth(2);
	gr->SetMarkerStyle(21);
	gr->SetMarkerSize(0.6);
	
	return gr;
	
}

TGraphAsymmErrors *plotEffPt(RooDataSet *a, int aa) {
    const RooArgSet *set = a->get();
    RooRealVar *xAx = (RooRealVar*)set->find("pt");
    RooRealVar *eff = (RooRealVar*)set->find("efficiency");

    const int nbins = xAx->getBinning().numBins();

    double tx[nbins], txhi[nbins], txlo[nbins];
    double ty[nbins], tyhi[nbins], tylo[nbins];

    for (int i=0; i<nbins; i++) {
        a->get(i);
        ty[i] = eff->getVal();
        tx[i] = xAx->getVal();
        txhi[i] = fabs(xAx->getErrorHi());
        txlo[i] = fabs(xAx->getErrorLo());
        tyhi[i] = fabs(eff->getErrorHi());
        tylo[i] = fabs(eff->getErrorLo());
    }

    cout<<"NBins : "<<nbins<<endl;

    const double *x = tx;
    const double *xhi = txhi;
    const double *xlo = txlo;
    const double *y = ty;
    const double *yhi = tyhi;
    const double *ylo = tylo;

    TGraphAsymmErrors *b = new TGraphAsymmErrors();
    if(aa == 1) {
        *b = TGraphAsymmErrors(nbins,x,y,xlo,xhi,ylo,yhi);
    }
    if(aa == 0) {
        *b = TGraphAsymmErrors(nbins,x,y,0,0,ylo,yhi);
    }
    b->SetMaximum(1.1);
    b->SetMinimum(0.0);
    b->SetMarkerStyle(20);
    b->SetMarkerColor(kRed+2);
    b->SetMarkerSize(1.0);
    b->SetTitle("");
    b->GetXaxis()->SetTitleSize(0.05);
    b->GetYaxis()->SetTitleSize(0.05);
    b->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    b->GetYaxis()->SetTitle("Efficiency");
    b->GetXaxis()->CenterTitle();
    //b->Draw("apz");

    for (int i=0; i<nbins; i++) {
        cout << x[i] << " " << y[i] << " " << yhi[i] << " " << ylo[i] << endl;
    }

    return b;

}

void FitterUtils::initiateParams(int nGenSignalZeroGamma, int nGenSignalOneGamma, int nGenSignalTwoGamma, RooRealVar const& expoConstGen, RooRealVar& nSignal, RooRealVar& nPartReco, 
      RooRealVar& nComb, RooRealVar& fracZero, RooRealVar& fracOne, RooRealVar& expoConst, RooRealVar&  nJpsiLeak, bool constPartReco, RooRealVar const& fracPartRecoSigma)
{
   TRandom rand;
   rand.SetSeed();

   int nGenSignal = nGenSignalZeroGamma + nGenSignalOneGamma + nGenSignalTwoGamma;

   double nGenSignal2;
   double nGenPartReco2;
   if(!constPartReco)
   {
      nGenSignal2 = rand.Uniform(nGenSignal-5*sqrt(nGenSignal), nGenSignal+5*sqrt(nGenSignal));
      nGenPartReco2 = rand.Uniform(nGenPartReco-5*sqrt(nGenPartReco), nGenPartReco+5*sqrt(nGenPartReco));
   }
   if(constPartReco)
   { 
      double nGenSigPartReco( nGenSignal+nGenPartReco );
      double nGenSigPartReco2( rand.Uniform( nGenSigPartReco-5*sqrt(nGenSigPartReco), nGenSigPartReco+5*sqrt(nGenSigPartReco) ) );
      double fracPartReco1( nGenPartReco/(1.*nGenSignal));
      double fracPartReco2( rand.Uniform(fracPartReco1-5*fracPartRecoSigma.getVal(), fracPartReco1+5*fracPartRecoSigma.getVal()) ); 

      nGenPartReco2 = fracPartReco2*nGenSigPartReco2 / (1+fracPartReco2); 
      nGenSignal2 = nGenSigPartReco2 / (1+fracPartReco2); 
   }
   double nGenComb2 = rand.Uniform(nGenComb-5*sqrt(nGenComb), nGenComb+5*sqrt(nGenComb));
   double nGenJpsiLeak2 = rand.Uniform(nGenJpsiLeak-5*sqrt(nGenJpsiLeak), nGenJpsiLeak+5*sqrt(nGenJpsiLeak));


   nSignal.setVal(nGenSignal2);
   nSignal.setRange(TMath::Max(0.,nGenSignal2-10.*sqrt(nGenSignal)) , nGenSignal2+10*sqrt(nGenSignal));

   nPartReco.setVal(nGenPartReco2);
   nPartReco.setRange(TMath::Max(0.,nGenPartReco2-10.*sqrt(nGenPartReco)), nGenPartReco2+10*sqrt(nGenPartReco));


   nComb.setVal(nGenComb2);
   nComb.setRange(TMath::Max(0.,nGenComb2-10.*sqrt(nGenComb)), nGenComb2+10*sqrt(nGenComb));

   nJpsiLeak.setVal(nGenJpsiLeak2);
   nJpsiLeak.setRange(TMath::Max(0., nGenJpsiLeak2-10*sqrt(nGenJpsiLeak)), nGenJpsiLeak2+10*sqrt(nGenJpsiLeak));

   double fracGenZero(nGenSignalZeroGamma/(1.*nGenSignal));
   double fracGenOne(nGenSignalOneGamma/(1.*nGenSignal));

   fracZero.setVal(rand.Gaus(fracGenZero, sqrt(nGenSignalZeroGamma)/(1.*nGenSignal))) ;
   fracZero.setRange(0., 1.);
   fracOne.setVal(rand.Gaus(fracGenOne, sqrt(nGenSignalOneGamma)/(1.*nGenSignal))) ;
   fracOne.setRange(0., 1.);

   expoConst.setVal(rand.Uniform( expoConstGen.getVal() - 5*expoConstGen.getError(), expoConstGen.getVal() + 5*expoConstGen.getError() ) );
   expoConst.setRange( expoConstGen.getVal() - 10*expoConstGen.getError(), expoConstGen.getVal() + 10*expoConstGen.getError() );
}

MakeBiasStudy::MakeBiasStudy() {
  int Nmodels = 8;
  //RooRealVar* mass = ws->var("mass");
  RooRealVar *mass = new RooRealVar("mass","mass", 100,180);
  RooRealVar *nBkgTruth = new RooRealVar("TruthNBkg","", 0,1e9);
  //  RooAbsData* realData = ws->data("Data_Combined")->reduce( Form("evtcat==evtcat::%s",cat.Data()) );

  double Bias[Nmodels][Nmodels];
  double BiasE[Nmodels][Nmodels];
  MakeAICFits MakeAIC_Fits;
  for(int truthType = 0; truthType < Nmodels; truthType++){

    RooAbsPdf *truthPdf = MakeAIC_Fits.getBackgroundPdf(truthType,mass);
    RooExtendPdf *truthExtendedPdf = new RooExtendPdf("truthExtendedPdf","",*truthPdf,*nBkgTruth);
    //truthExtendedPdf.fitTo(*realData,RooFit::Strategy(0),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));
    //truthExtendedPdf.fitTo(*realData,RooFit::Strategy(2),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));
    
    double BiasWindow = 2.00;
    mass->setRange("biasRegion", mh-BiasWindow, mh+BiasWindow);
    double TruthFrac = truthExtendedPdf->createIntegral(mass,RooFit::Range("biasRegion"),RooFit::NormSet(*mass))->getVal();
    double NTruth = TruthFrac * nBkgTruth->getVal();
    double NTruthE = TruthFrac * nBkgTruth->getError();
    
    RooDataSet* truthbkg = truthPdf->generate(RooArgSet(*mass),nBkgTruth);

    for(int modelType = 0; modelType < Nmodels; modelType++){
      RooAbsPdf* ModelShape = MakeAIC_Fits.getBackgroundPdf(modelType,mass);
      RooRealVar *nBkgFit = new RooRealVar("FitNBkg", "", 0, 1e9);
      RooExtendPdf ModelExtendedPdf = new RooExtendPdf("ModelExtendedPdf", "",*ModelShape, *nBkgFit);
      ModelExtendedPdf.fitTo(truthbkg, RooFit::Strategy(0),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));
      ModelExtendedPdf.fitTo(truthbkg, RooFit::Strategy(2),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));

      double FitFrac = ModelExtendedPdf.createIntegral(mass,RooFit::Range("biasRegion"),RooFit::NormSet(mass))->getVal();
      double NFit = FitFrac * nBkgFit->getVal();
      double NFitE = FitFrac * nBkgFit->getError();

      Bias[truthType][modelType] = fabs(NFit - NTruth);
      BiasE[truthType][modelType] = fabs(NFitE - NTruthE);
    }

    
  }
  
  for(int i = 0; i < Nmodels; i++) { 
    std::cout <<  "===== Truth Model : " << MakeBiasStudy::Category(i) << " ===== " << std::endl; 
    for (int j = 0; j < Nmodels; j++) { 
      std::cout << "Fit Model: " << MakeBiasStudy::Category(j) << "  , Bias = " << Bias[i][j] << " +/- " << BiasE[i][j] << std::endl; 
    }   
  }
  

}

void MakeSpinPlots::MakeChannelComp(const char* mcType){
  TGraphErrors graph(catNames.size());

  RooRealVar mu("mu","",-50,50);

  float exp = ws->data(Form("%s_Combined",mcType))->sumEntries();///total * 607*lumi/12.;

  TH1F frame("frame","",catNames.size(),0,catNames.size());

  //graph.GetXaxis()->SetNdivisions(catNames.size());

  float min=99999,max=-99999;
  for(int i=0;i<catNames.size();i++){
    RooRealVar *ind = ws->var( Form("Data_%s_INDFIT_%s_Nsig",mcType, catNames.at(i).Data()) );
    RooRealVar *f = ws->var( Form("Data_%s_FULLFIT_%s_fsig",mcType, catNames.at(i).Data()) );
    float mu = ind->getVal()/exp/f->getVal();
    float muE = ind->getError()/exp/f->getVal();
    graph.SetPoint(i,i+0.5,mu);
    graph.SetPointError(i,0,muE);

    if(mu-muE < min) min = mu-muE;
    if(mu+muE > max) max = mu+muE;
    
    //graph.GetXaxis()->SetBinLabel(i+1,catNames.at(i));
  }

  TF1 fit("fit","[0]",0+0.5,catNames.size()+0.5);
  graph.Fit(&fit,"MNE");
  float mean  = fit.GetParameter(0);
  float meanE = fit.GetParError(0);
  
  frame.SetAxisRange(min-1.5,max+0.5,"Y");
  frame.SetYTitle("Fitted #sigma/#sigma_{SM}");
  frame.SetXTitle("Category");
  TCanvas cv;
  frame.Draw();

  TBox err(0,mean-meanE,catNames.size(),mean+meanE);
  err.SetFillColor(kGreen);

  frame.Draw();
  err.Draw("SAME");
  TLine mLine(0,mean,catNames.size(),mean);
  mLine.Draw("SAME");
  graph.Draw("PSAME");
  

  cv.SaveAs(basePath+Form("/ChannelComp_%s_%s.png",mcType,outputTag.Data()));
  cv.SaveAs(basePath+Form("/ChannelComp_%s_%s.pdf",mcType,outputTag.Data()));
}

double getSigmaError(RooWorkspace *w) 
{
  using namespace RooFit;
  // Access saved Sigma CB values and error from WorkSpace
  RooRealVar *mySavedSigma = w->var("sigma_{CB}");
  RooRealVar *mySavedDeltaM = w->var("#Deltam_{CB}");
  Double_t sigma = mySavedSigma->getVal();
  Double_t sigmaerror = mySavedSigma->getError();
  Double_t DeltaM = mySavedDeltaM->getVal();
  Double_t DeltaMerror = mySavedDeltaM->getError();
  Double_t Sigpercent = 100 * sigma / (91.187 + DeltaM);
  Double_t Sigpercenterror = 100 * pow(pow(sigmaerror / (91.198 + DeltaM), 2) + pow((sigma * DeltaMerror) / pow(DeltaM + 91.198, 2), 2), .5);
  return Sigpercenterror;
} // plotFitOn(..)

void testResolution() { 
  Prep();
  RooRealVar    lXVar  ("XVar","mass(GeV/c^{2})",100,60,150); lXVar.setBins(1000);
  RooRealVar l1Sigma("sigma1","sigma1",1.6 ,0.,15.);  //l1Sigma.setConstant(kTRUE);
  RooRealVar l2Sigma("sigma2","sigma2",1.6,0.,15.);   l2Sigma.setConstant(kTRUE);
  RooRealVar l3Sigma("sigma3","sigma3",2.90,0.,35.);  l3Sigma.setConstant(kTRUE);
  RooRealVar lN     ("n"     ,"n"     ,1.00,-15,15.); lN.setConstant(kTRUE);
  RooRealVar lExp   ("exp"   ,"exp"   ,-0.003,-15,15.); //lExp.setConstant(kTRUE);
  RooRealVar lR0Mean("xmean","xmean",0,-10,10);    lR0Mean.setConstant(kTRUE);
  RooRealVar lR1Mean("mean","mean",90.8,60,150);   //lR1Mean.setConstant(kTRUE);
  RooVoigtianShape     lGAdd("Add","Add",lXVar,lR1Mean,l1Sigma,l2Sigma,lN,l3Sigma,true);
  
  RooRealVar    lSPar  ("SPar","SPar", 1.,0., 2.);
  RooFormulaVar lXShift("uparshift","@0*@1",RooArgList(lXVar,lSPar));
  TH1F *lMass = getMass(0,-5,5,-1.5,1.5);
  RooDataHist *lMHist = new RooDataHist("M" ,"M" ,RooArgSet(lXVar),lMass);
  RooHistPdf  *lMPdf  = new RooHistPdf ("MH","MH",lXShift,lXVar,*lMHist,5); 
  RooRealVar lGSigma("gsigma","gsigma",1.6 ,0.,15.);  
  RooGaussian   lGaus1("gaus1","gaus1",lXVar,lR0Mean,lGSigma);
  RooFFTConvPdf  lConv("Conv","Conv",lXVar,*lMPdf,lGaus1);

  RooDataSet *lData = new RooDataSet("crap","crap",RooArgSet(lXVar)); 
  fillData(lData,lXVar,0,-5,5,-1.5,1.5);
  lConv.fitTo(*lData,Strategy(2));
  lGSigma.setVal(lGSigma.getVal());
  lSPar.setVal(lSPar.getVal()*1.01);
  //cout << "=====> Check " << l1Sigma.getVal() << " --- " << lR1Mean.getVal() << "----" << lSPar.getVal() << endl;
  PEs(&lConv,&lGAdd,2000,500,lXVar,l1Sigma,lR1Mean,lSPar,lGSigma);
  lData = 0;
  Plot(&lConv,&lGAdd,2000,50000,lXVar,l1Sigma,lR1Mean,lSPar,lGSigma,lData);
}

void simplePrintResults() {
    vector<string> filenames;
    filenames.push_back("Output/Test2/result/FIT_DATA_Psi2SJpsi_PPPrompt_Bkg_SecondOrderChebychev_pt65300_rap016_cent0200_262620_263757.root");

    const char* parname = "N_Jpsi_PP";

    vector<string>::iterator it = filenames.begin();
    for (it; it<filenames.end(); it++) {
        TFile *f = new TFile(it->c_str());
        if (!f) {
            cout << "Error, " << *it << " not found" << endl;
            continue;
        }
        RooWorkspace *ws = (RooWorkspace*) f->Get("workspace");
        if (!ws) {
            cout << "Error, workspace not found in " << *it << endl;
            continue;
        }

        RooRealVar *var = ws->var(parname);
        if (!ws) {
            cout << "Error, variable " << parname << " not found in " << *it << endl;
            continue;
        }
        cout << *it << " " << var->getVal() << " +- " << var->getError() << endl;
    }
}

void PDF_GLWADS_DKDpi_K3pi::setObservables(config c)
{
	switch(c)
	{
		case truth:{
					   setObservablesTruth();
					   break;
				   }
		case toy:{
					 setObservablesToy();
					 break;
				 }
		case lumi1fb:{
						 obsValSource = "1fb-1, ExpNll/sept2012K3PIResult.root";
						 TString File = this->dir+"/ExpNll/sept2012K3PIResult.root";
						 TFile *fr = TFile::Open(File);
						 RooFitResult *r = (RooFitResult*)fr->Get("fitresult_model_reducedData_binned");
						 assert(r);
						 TIterator* it = observables->createIterator();
						 while ( RooRealVar* pObs = (RooRealVar*)it->Next() )
						 {
							 RooRealVar* pRes = (RooRealVar*)r->floatParsFinal().find(obsTmkToMalcolm(pObs->GetName()));
							 pObs->setVal(pRes->getVal());
						 }
						 fr->Close();
						 delete r;
						 delete fr;
						 break;
					 }
		case lumi3fb:{
						 obsValSource = "3fb-1 ANA v7 unblind"; // https://twiki.cern.ch/twiki/pub/LHCbPhysics/B2D0K/LHCb-ANA-2014-071-v7.pdf (see Vavas email 04/08/15)

						 // these get transformed over from the new inputs using ExpNll/transportGLWADS_new_to_old.py
						 // in the case of the DK only (robust) combination some of the observables don't exist
						 // usemap as the temp store
						 std::map< TString, double > vals;
						 vals["rkp_k3pi_obs"]      =  0.0793;
						 vals["afav_dk_k3pi_obs"]  =  -0.0004;
						 vals["afav_dpi_k3pi_obs"] =  0.0;
						 vals["rp_dk_k3pi_obs"]    =  0.018369;
						 vals["rm_dk_k3pi_obs"]    =  0.009611;
						 vals["rp_dpi_k3pi_obs"]   =  0.003683;
						 vals["rm_dpi_k3pi_obs"]   =  0.003857;

						 // now can loop the observables and set the values
						 TIterator* it = observables->createIterator();
						 while ( RooRealVar* pObs = (RooRealVar*)it->Next() ){
							 pObs->setVal(vals[pObs->GetName()]);
						 }

						 vals.clear();
						 break;
				}
		default:{
					cout << "PDF_GLWADS_DKDpi_K3pi::setObservables() : ERROR : config "+ConfigToTString(c)+" not found." << endl;
					exit(1);
				}
	}
}

///
/// Set each parameter in setMe to the value found in values.
/// Do nothing if parameter is not found in values.
///
void Utils::setParameters(const RooAbsCollection* setMe, const RooAbsCollection* values)
{
	TIterator* it = setMe->createIterator();
	while ( RooRealVar* p = (RooRealVar*)it->Next() ){
		RooRealVar *var = (RooRealVar*)values->find(p->GetName());
		if ( var ) p->setVal(var->getVal());
	}
	delete it;
}

// grab the initial parameters and errors for making pull distributions:
// Take these from a fit file to the data themselves 
void fillInitialParams(RooArgSet *args, std::map<std::string, std::pair<double,double> > &vals){
	
	 TIterator* iter(args->createIterator());
         for (TObject *a = iter->Next(); a != 0; a = iter->Next()) {
                 RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);      
                 std::string name = rrv->GetName();
		 std::pair<double,double> valE(rrv->getVal(),rrv->getError());
		 vals.insert( std::pair<std::string,std::pair<double ,double> > (name,valE)) ;
	 }
	
}

void GetNominalValueNuisancePara(){
    TIterator *it = mc->GetNuisanceParameters()->createIterator();
    RooRealVar *var = NULL;
    if (MapNuisanceParamNom.size() > 0) MapNuisanceParamNom.clear();
    std::cout << "Nuisance parameter names and values" << std::endl;
    while ((var = (RooRealVar*)it->Next()) != NULL){
      const double val = var->getVal();
      MapNuisanceParamNom[(string)var->GetName()] = val;
    }
    return;
 }

TH2F *plotEff2D(RooDataSet *a, TString b) {
    const RooArgSet *set = a->get();
    RooRealVar *yAx = (RooRealVar*)set->find("pt");
    RooRealVar *xAx = (RooRealVar*)set->find("eta");
    RooRealVar *eff = (RooRealVar*)set->find("efficiency");

//    const int xnbins = xAx->getBinning().numBins();
//    const int ynbins = yAx->getBinning().numBins();

    //double xbins[] = {-2.4, -1.6, -0.8, 0.0, 0.8, 1.6, 2.4};
    //double ybins[] = {0, 2, 3, 5, 8, 10, 20};

    const double *xvbins = xAx->getBinning().array();
    const double *yvbins = yAx->getBinning().array();

    TH2F* h = new TH2F(b, "", xAx->getBinning().numBins(), xvbins, yAx->getBinning().numBins(), yvbins);

    gStyle->SetPaintTextFormat("5.2f");
    gStyle->SetPadRightMargin(0.12);
    gStyle->SetPalette(1);
    h->SetOption("colztexte");
    h->GetZaxis()->SetRangeUser(-0.001,1.001);
    h->SetStats(kFALSE);
    h->GetYaxis()->SetTitle("p_{T} [GeV/c]");
    h->GetXaxis()->SetTitle("#eta");
    h->GetXaxis()->CenterTitle();
    h->GetYaxis()->CenterTitle();
    h->GetXaxis()->SetTitleSize(0.05);
    h->GetYaxis()->SetTitleSize(0.05);
    h->GetYaxis()->SetTitleOffset(0.8);
    h->GetXaxis()->SetTitleOffset(0.9);
    for(int i=0; i<a->numEntries(); i++) {
        a->get(i);
        h->SetBinContent(h->FindBin(xAx->getVal(), yAx->getVal()), eff->getVal());
        h->SetBinError(h->FindBin(xAx->getVal(), yAx->getVal()), (eff->getErrorHi()-eff->getErrorLo())/2.);
    }

    return h;

}