C++ RooFitResult::floatParsFinal方法代码示例

MCMCInterval * Tprime::GetMcmcInterval(ModelConfig mc,
                                       double conf_level,
                                       int n_iter,
                                       int n_burn,
                                       double left_side_tail_fraction,
                                       int n_bins) {
    //
    // Bayesian MCMC calculation using arbitrary ModelConfig
    // Want an efficient proposal function, so derive it from covariance
    // matrix of fit
    //

    RooAbsData * _data = data;
    //RooAbsData * _data = pWs->data("obsData");
    //RooStats::ModelConfig * _mc = (RooStats::ModelConfig *)pWs->genobj("ModelConfig");
    RooStats::ModelConfig * _mc = GetModelConfig();
    _mc->Print();

    //RooFitResult * fit = pWs->pdf("model_tprime")->fitTo(*_data,Save());
    RooFitResult * fit = _mc->GetPdf()->fitTo(*_data,Save());
    ProposalHelper ph;
    ph.SetVariables((RooArgSet&)fit->floatParsFinal());
    ph.SetCovMatrix(fit->covarianceMatrix());
    ph.SetUpdateProposalParameters(kTRUE); // auto-create mean vars and add mappings
    ph.SetCacheSize(100);
    ProposalFunction * pf = ph.GetProposalFunction();

    //delete pf;
    //pf = new SequentialProposal();

    MCMCCalculator mcmc( *_data, mc );
    mcmc.SetConfidenceLevel(conf_level);
    mcmc.SetNumIters(n_iter);          // Metropolis-Hastings algorithm iterations
    mcmc.SetProposalFunction(*pf);
    mcmc.SetNumBurnInSteps(n_burn); // first N steps to be ignored as burn-in
    mcmc.SetLeftSideTailFraction(left_side_tail_fraction);
    mcmc.SetNumBins(n_bins);

    //mcInt = mcmc.GetInterval();
    try {
        mcInt = mcmc.GetInterval();
    } catch ( std::length_error &ex) {
        mcInt = 0;
    }

    //std::cout << "!!!!!!!!!!!!!! interval" << std::endl;
    if (mcInt == 0) std::cout << "No interval found!" << std::endl;

    delete fit;
    delete pf;

    return mcInt;
}

  void   printResult(){      
    if (fitresult!=NULL){ 
      printf("%s\n",
	     "==================================================+++" );
      fitresult->Print("v")  ;   
      fitresult->floatParsFinal().Print("s") ;
      printf("%s\nfit range: %5.1f %5.1f\n",
	     "==================================================+++",
               min, max );
      //RooRealVar* par1_fitresult = (RooRealVar*) fitresult->floatParsFinal()->find("par1") 
      //par1_fitresult->GetAsymErrorHi() ; // etc... 
    }
  }//printResult

MCMCInterval * TwoBody::GetMcmcInterval_OldWay(ModelConfig mc,
					double conf_level,
					int n_iter,
					int n_burn,
					double left_side_tail_fraction,
					int n_bins){
  // use MCMCCalculator  (takes about 1 min)
  // Want an efficient proposal function, so derive it from covariance
  // matrix of fit
  
  RooFitResult* fit = ws->pdf("model")->fitTo(*data,Save());
  ProposalHelper ph;
  ph.SetVariables((RooArgSet&)fit->floatParsFinal());
  ph.SetCovMatrix(fit->covarianceMatrix());
  ph.SetUpdateProposalParameters(kTRUE); // auto-create mean vars and add mappings
  ph.SetCacheSize(100);
  ProposalFunction* pf = ph.GetProposalFunction();
  
  MCMCCalculator mcmc( *data, mc );
  mcmc.SetConfidenceLevel(conf_level);
  mcmc.SetNumIters(n_iter);          // Metropolis-Hastings algorithm iterations
  mcmc.SetProposalFunction(*pf);
  mcmc.SetNumBurnInSteps(n_burn); // first N steps to be ignored as burn-in
  mcmc.SetLeftSideTailFraction(left_side_tail_fraction);
  mcmc.SetNumBins(n_bins);
  
//mcInt = mcmc.GetInterval();
  try {
    mcInt = mcmc.GetInterval();
  } catch ( std::length_error &ex) {
    mcInt = 0;
  }

  //std::cout << "!!!!!!!!!!!!!! interval" << std::endl;
  if (mcInt == 0) std::cout << "No interval found!" << std::endl;
  
  return mcInt;
}

///
/// Fit a pdf to the minimum, but keep angular parameters in a range of
/// [0,2pi]. If after an initial fit, a parameter has walked outside this
/// interval, add multiples of 2pi to bring it back. Then, refit.
/// All variables that have unit 'rad' are taken to be angles.
///
RooFitResult* Utils::fitToMinBringBackAngles(RooAbsPdf *pdf, bool thorough, int printLevel)
{
	countAllFitBringBackAngle++;
	RooFitResult* r = fitToMin(pdf, thorough, printLevel);
	bool refit = false;
	TIterator* it = r->floatParsFinal().createIterator();
	while ( RooRealVar* p = (RooRealVar*)it->Next() ){
		if ( ! isAngle(p) ) continue;
		if ( p->getVal()<0.0 || p->getVal()>2.*TMath::Pi() ){
			RooArgSet *pdfPars = pdf->getParameters(RooArgSet());
			RooRealVar *pdfPar = (RooRealVar*)pdfPars->find(p->GetName());
			pdfPar->setVal(bringBackAngle(p->getVal()));
			refit = true;
			delete pdfPars;
		}
	}
	if ( refit ){
		countFitBringBackAngle++;
		delete r;
		r = fitToMin(pdf, thorough, printLevel);
	}
	delete it;
	return r;
}

void Fit3D::plotFitProjection(
    const RooRealVar &independant_variable,
    const RooDataSet &data,
    const RooFitResult& fit,
    const RooAbsPdf &model,
    const RooAbsPdf &bs_pdf,
    const RooAbsPdf &bd_pdf,
    const RooAbsPdf &cw_pdf,
    const RooAbsPdf &ww_pdf,
    const RooAbsPdf &cn_pdf,
    const TString &filename)
{
  RooPlot* frame = independant_variable.frame();
  TString frame_title = "Fit Projection on ";
  frame_title.Append(independant_variable.GetTitle());
  frame->SetTitle(frame_title);
  
  RooRealVar* bs_fit = (RooRealVar*) fit.floatParsFinal().find("n_bs_pp");
  RooRealVar* bd_fit = (RooRealVar*) fit.floatParsFinal().find("n_bd_pp");
  RooRealVar* cw_fit = (RooRealVar*) fit.floatParsFinal().find("n_cw_pp");
  RooRealVar* ww_fit = (RooRealVar*) fit.floatParsFinal().find("n_ww_pp");
  RooRealVar* cn_fit = (RooRealVar*) fit.floatParsFinal().find("n_cn_pp");
  
  if (!bs_fit) {
    bs_fit = (RooRealVar*) fit.floatParsFinal().find("n_bs_nn");
    bd_fit = (RooRealVar*) fit.floatParsFinal().find("n_bd_nn");
    cw_fit = (RooRealVar*) fit.floatParsFinal().find("n_cw_nn");
    ww_fit = (RooRealVar*) fit.floatParsFinal().find("n_ww_nn");
    cn_fit = (RooRealVar*) fit.floatParsFinal().find("n_cn_nn");
  }
  if (!bs_fit) {
    // Error. Quit while ahead.
    cout << "Error in plotFitAccuracy(): "
         << "Cannot find fit variables. Check names are valid."
         << endl;
    return;
  }
  
  data.plotOn(frame, RooFit::Name("data"));
  // model.plotOn(frame, RooFit::Name("model"), RooFit::LineColor(kBlue));
  bs_pdf.plotOn(frame,
      RooFit::Normalization(bs_fit->getVal(), RooAbsReal::NumEvent),
      RooFit::LineStyle(kDashed),
      RooFit::LineWidth(1),
      RooFit::LineColor(kYellow + 2));
  data.plotOn(frame, RooFit::Cut(bs_events_cut_),
      RooFit::LineColor(kYellow),
      RooFit::MarkerStyle(kFullDotMedium));
  
  bd_pdf.plotOn(frame,
      RooFit::Normalization(bd_fit->getVal(), RooAbsReal::NumEvent),
      RooFit::LineStyle(kDashed),
      RooFit::LineWidth(1),
      RooFit::LineColor(kRed + 2));
  data.plotOn(frame, RooFit::Cut(bd_events_cut_),
      RooFit::LineColor(kRed),
      RooFit::MarkerStyle(kFullDotMedium));
  cw_pdf.plotOn(frame,
      RooFit::Normalization(cw_fit->getVal(), RooAbsReal::NumEvent),
      RooFit::LineStyle(kDashed),
      RooFit::LineWidth(1),
      RooFit::LineColor(kGreen + 2));
  data.plotOn(frame, RooFit::Cut(cw_events_cut_),
      RooFit::LineColor(kGreen),
      RooFit::MarkerStyle(kFullDotMedium));
  ww_pdf.plotOn(frame,
      RooFit::Normalization(ww_fit->getVal(), RooAbsReal::NumEvent),
      RooFit::LineStyle(kDashed),
      RooFit::LineWidth(1),
      RooFit::LineColor(kBlue + 2));
  data.plotOn(frame, RooFit::Cut(ww_events_cut_),
      RooFit::LineColor(kBlue),
      RooFit::MarkerStyle(kFullDotMedium));
  cn_pdf.plotOn(frame,
      RooFit::Normalization(cn_fit->getVal(), RooAbsReal::NumEvent),
      RooFit::LineStyle(kDashed),
      RooFit::LineWidth(1),
      RooFit::LineColor(kCyan + 2));
  data.plotOn(frame, RooFit::Cut(cn_events_cut_),
      RooFit::LineColor(kCyan),
      RooFit::MarkerStyle(kFullDotMedium));
      
  TCanvas* c1 = new TCanvas("c1", "Projection", 200, 10, 700, 500);
  frame->Draw();
  c1->Print(output_path_ + filename);
}

//.........这里部分代码省略.........
  pdf_KD_sm.plotOn(kdframe1,LineColor(kBlue));
  model.plotOn(kdframe1,LineColor(kRed));
  kdframe1->Draw();
  c->SaveAs("model1DPlot.png");
  c->SaveAs("model1DPlot.eps");

  double fa3Val=-99;
  if (mySample == kScalar_fa3_0)
    fa3Val=0.;
  else if (mySample == kScalar_fa3_1)
    fa3Val=1;
  else{
    cout<<"fa3Val not correct!"<<endl;
      return 0;
  }

  TCanvas* c = new TCanvas("modelPlot","modelPlot",400,400);
  rdh_KD_ps.plotOn(kdframe1,LineColor(kBlack));
  pdf_KD_ps.plotOn(kdframe1,LineColor(kBlack));
  rdh_KD_sm.plotOn(kdframe1,LineColor(kBlue));
  pdf_KD_sm.plotOn(kdframe1,LineColor(kBlue));
  model.plotOn(kdframe1,LineColor(kRed));
  kdframe1->Draw();

  
  TChain* myChain = new TChain("newTree");
  myChain->Add(inputFileNames[mySample]);
  
  if(!myChain || myChain->GetEntries()<=0) {
    cout<<"error in the tree"<<endl;
    return 0;
  }
  
   RooDataSet* data = new RooDataSet("data","data",myChain,RooArgSet(*kd),"");

    cout << "Number of events in data: " << data->numEntries() << endl;
  
  // initialize tree to save toys to 
  TTree* results = new TTree("results","toy results");
  
  double fa3,fa3Error, fa3Pull;
  double fa2,fa2Error, fa2Pull;
  double phia2, phia2Error, phia2Pull;
  double phia3, phia3Error, phia3Pull;

  results->Branch("fa3",&fa3,"fa3/D");
  results->Branch("fa3Error",&fa3Error,"fa3Error/D");
  results->Branch("fa3Pull",&fa3Pull,"fa3Pull/D");

  //---------------------------------

  RooDataSet* toyData;
  int embedTracker=0;
  RooArgSet *tempEvent;

  RooFitResult *toyfitresults;
  RooRealVar *r_fa3;

  for(int i = 0 ; i<nToys ; i++){
    cout <<i<<"<-----------------------------"<<endl;
    //if(toyData) delete toyData;
    toyData = new RooDataSet("toyData","toyData",RooArgSet(*kd));

    if(nEvts+embedTracker > data->sumEntries()){
      cout << "Playground::generate() - ERROR!!! Playground::data does not have enough events to fill toy!!!!  bye :) " << endl;
      toyData = NULL;
      abort();
      return 0;
    }

    for(int iEvent=0; iEvent<nEvts; iEvent++){
      if(iEvent==1)
	cout << "generating event: " << iEvent << " embedTracker: " << embedTracker << endl;
      tempEvent = (RooArgSet*) data->get(embedTracker);
      toyData->add(*tempEvent);
      embedTracker++;
    }

    toyfitresults =model.fitTo(*toyData,Save());
    cout<<toyfitresults<<endl;
    r_fa3 = (RooRealVar *) toyfitresults->floatParsFinal().find("fa3");

    fa3 = r_fa3->getVal();
    fa3Error = r_fa3->getError();
    fa3Pull = (r_fa3->getVal() - fa3Val) / r_fa3->getError();

    // fill TTree
    results->Fill();
    //model.clear();
  }

  char nEvtsString[100];
  sprintf(nEvtsString,"_%iEvts",nEvts);

  // write tree to output file (ouputFileName set at top)
  TFile *outputFile = new TFile("embeddedToysVBF_fa3Corr_"+sampleName[mySample]+nEvtsString+".root","RECREATE");
  results->Write();
  outputFile->Close();

}

void plotResolution() {

  TStyle *mystyle = RooHZZStyle("ZZ");
  mystyle->cd();

  double binedgesZ[5] = {20,30,40,50,70};
  double binedgesJPsi[4] = {7,10,20,40};
  TGraphAsymmErrors gScaleZ[2][2];
  TGraphAsymmErrors gResoZ[2][2];
  TGraphAsymmErrors gScaleJPsi[2];
  TGraphAsymmErrors gResoJPsi[2];

  // Z->ee
  for(int ipt=0; ipt<4; ++ipt) {
    for(int ieta=0; ieta<2; ++ieta) {
      for(int ir9=0; ir9<2; ++ir9) {

	cout << "Analyzing pt bin = " << ipt << ", eta bin: " << ieta << "  and r9 bin: " << ir9 << endl;

	stringstream mcfile, datafile;
	mcfile << "mcZ2012_PtBin" << ipt << "_EtaBin" << ieta << "_R9Bin" << ir9 << ".root";
	datafile << "dataZ2012_PtBin" << ipt << "_EtaBin" << ieta << "_R9Bin" << ir9 << ".root";
    
	TFile *tmcfile = TFile::Open(mcfile.str().c_str());
	RooFitResult *mcfr = (RooFitResult*)tmcfile->Get("fitres");
	float mcDM = ((RooRealVar*)(mcfr->floatParsFinal().find("#Deltam_{CB}")))->getVal();
	float mcDM_err = ((RooRealVar*)(mcfr->floatParsFinal().find("#Deltam_{CB}")))->getError();
	float mcS = ((RooRealVar*)(mcfr->floatParsFinal().find("#sigma_{CB}")))->getVal();
	float mcS_err = ((RooRealVar*)(mcfr->floatParsFinal().find("#sigma_{CB}")))->getError();
	mcS_err = effRmsFromSigmaCB(mcDM,mcDM_err,mcS,mcS_err);
	

	TFile *tdatafile = TFile::Open(datafile.str().c_str());
	RooFitResult *datafr = (RooFitResult*)tdatafile->Get("fitres");
	float dataDM = ((RooRealVar*)(datafr->floatParsFinal().find("#Deltam_{CB}")))->getVal();
	float dataDM_err = ((RooRealVar*)(datafr->floatParsFinal().find("#Deltam_{CB}")))->getError();
	float dataS = ((RooRealVar*)(datafr->floatParsFinal().find("#sigma_{CB}")))->getVal();
	float dataS_err = ((RooRealVar*)(datafr->floatParsFinal().find("#sigma_{CB}")))->getError();
	dataS_err = effRmsFromSigmaCB(dataDM,dataDM_err,dataS,dataS_err);
	
	float rM = (dataDM-mcDM)/(dataDM + 91.19);
	float rM_err = rM * sqrt(dataDM_err*dataDM_err + mcDM_err*mcDM_err);
	float delta = dataS-mcS;
	float delta_err = sqrt(dataS_err*dataS_err+mcS_err*mcS_err);
	float rS = delta/(mcS);
	cout << "rS = " << rS << endl;
	float rS_err = fabs(rS) * sqrt(pow(delta_err/delta,2) + pow(mcS_err/mcS,2));
	
	float bincenter=(binedgesZ[ipt+1]+binedgesZ[ipt])/2.;
	// add some offset not to overlap points
	bincenter += (ieta==0) ? 2. : -2.;
	bincenter += (ir9==0) ? 1. : -1.;
	
	float binerrup=binedgesZ[ipt+1]-bincenter;
	float binerrdn=bincenter-binedgesZ[ipt];

	gScaleZ[ieta][ir9].SetPoint(ipt,bincenter,rM);
	gScaleZ[ieta][ir9].SetPointError(ipt,binerrdn,binerrup,rM_err,rM_err);
	gResoZ[ieta][ir9].SetPoint(ipt,bincenter,rS);
	gResoZ[ieta][ir9].SetPointError(ipt,binerrdn,binerrup,rS_err,rS_err);
      }
    }
  }

  // JPsi->ee
  for(int ipt=0; ipt<3; ++ipt) {
    for(int ieta=0; ieta<1; ++ieta) {

	cout << "Analyzing pt bin = " << ipt << ", eta bin: " << ieta << endl;

	stringstream mcfile, datafile;
	mcfile << "mcJPsi2012_PtBin" << ipt << "_EtaBin" << ieta << ".root";
	datafile << "dataJPsi2012_PtBin" << ipt << "_EtaBin" << ieta << ".root";
    
	TFile *tmcfile = TFile::Open(mcfile.str().c_str());
	RooFitResult *mcfr = (RooFitResult*)tmcfile->Get("fitres");
	float mcDM = ((RooRealVar*)(mcfr->floatParsFinal().find("#Deltam_{CB}")))->getVal();
	float mcDM_err = ((RooRealVar*)(mcfr->floatParsFinal().find("#Deltam_{CB}")))->getError();
	float mcS = ((RooRealVar*)(mcfr->floatParsFinal().find("#sigma_{CB}")))->getVal();
	float mcS_err = ((RooRealVar*)(mcfr->floatParsFinal().find("#sigma_{CB}")))->getError();

	TFile *tdatafile = TFile::Open(datafile.str().c_str());
	RooFitResult *datafr = (RooFitResult*)tdatafile->Get("fitres");
	float dataDM = ((RooRealVar*)(datafr->floatParsFinal().find("#Deltam_{CB}")))->getVal();
	float dataDM_err = ((RooRealVar*)(datafr->floatParsFinal().find("#Deltam_{CB}")))->getError();
	float dataS = ((RooRealVar*)(datafr->floatParsFinal().find("#sigma_{CB}")))->getVal();
	float dataS_err = ((RooRealVar*)(datafr->floatParsFinal().find("#sigma_{CB}")))->getError();
	
	float rM = (dataDM-mcDM)/(dataDM + 3.096);
	float rM_err = rM * sqrt(dataDM_err*dataDM_err + mcDM_err*mcDM_err);
	float delta = dataS-mcS;
	float delta_err = sqrt(dataS_err*dataS_err+mcS_err*mcS_err);
	float rS = delta/(mcS);
	cout << "rS = " << rS << endl;
	float rS_err = fabs(rS) * sqrt(pow(delta_err/delta,2) + pow(mcS_err/mcS,2));
	
	float bincenter=(binedgesJPsi[ipt+1]+binedgesJPsi[ipt])/2.;
	// add some offset not to overlap points
	bincenter += (ieta==0) ? 0.5 : -0.5;
	
//.........这里部分代码省略.........

void Fit3D::plotFitAccuracy(
    const RooDataSet& mc_data,
    const RooFitResult& fit)
{
  fit.Print("v");

  double n_true_bs = mc_data.sumEntries("component == component::bs");
  double n_true_bd = mc_data.sumEntries("component == component::bd");
  double n_true_cw = mc_data.sumEntries("component == component::cw");
  double n_true_ww = mc_data.sumEntries("component == component::ww");
  double n_true_cn = mc_data.sumEntries("component == component::cn");

  RooRealVar* bs_fit = (RooRealVar*) fit.floatParsFinal().find("n_bs_pp");
  RooRealVar* bd_fit = (RooRealVar*) fit.floatParsFinal().find("n_bd_pp");
  RooRealVar* cw_fit = (RooRealVar*) fit.floatParsFinal().find("n_cw_pp");
  RooRealVar* ww_fit = (RooRealVar*) fit.floatParsFinal().find("n_ww_pp");
  RooRealVar* cn_fit = (RooRealVar*) fit.floatParsFinal().find("n_cn_pp");
  
  TString title("Fit Accuracy (N^{++}_{fit}-N^{++}_{true})");
  TString filename(output_path_ + "fit_accuracy_pp");
  if (!bs_fit) {
    bs_fit = (RooRealVar*) fit.floatParsFinal().find("n_bs_nn");
    bd_fit = (RooRealVar*) fit.floatParsFinal().find("n_bd_nn");
    cw_fit = (RooRealVar*) fit.floatParsFinal().find("n_cw_nn");
    ww_fit = (RooRealVar*) fit.floatParsFinal().find("n_ww_nn");
    cn_fit = (RooRealVar*) fit.floatParsFinal().find("n_cn_nn");
    title = TString("Fit Accuracy (N^{--}_{fit}-N^{--}_{true})");
    filename = TString(output_path_ + "fit_accuracy_nn");
  }
  if (!bs_fit) {
    // Error. Quit while ahead.
    cout << "Error in plotFitAccuracy(): "
         << "Cannot find fit variables. Check names are valid."
         << endl;
    return;
  }

  std::cout << n_true_bd << std::endl;
  std::cout << n_true_bs << std::endl;
  std::cout << n_true_cn << std::endl;
  std::cout << n_true_cw << std::endl;
  std::cout << n_true_ww << std::endl;
  
  TCanvas c1("c1", title, 200, 10, 700, 500);
  c1.SetGrid();
  double x[5] = {1, 2, 3, 4, 5};
  double y[5] = {
      bs_fit->getVal() - n_true_bs,
      bd_fit->getVal() - n_true_bd,
      cw_fit->getVal() - n_true_cw,
      ww_fit->getVal() - n_true_ww,
      cn_fit->getVal() - n_true_cn};
  double exl[5] = {0, 0, 0, 0, 0};
  double exh[5] = {0, 0, 0, 0, 0};
  double eyl[5] = {
      -bs_fit->getErrorLo(),
      -bd_fit->getErrorLo(),
      -cw_fit->getErrorLo(),
      -ww_fit->getErrorLo(),
      -cn_fit->getErrorLo()};
  double eyh[5] = {
      bs_fit->getErrorHi(),
      bd_fit->getErrorHi(),
      cw_fit->getErrorHi(),
      ww_fit->getErrorHi(),
      cn_fit->getErrorHi()};
  TGraphAsymmErrors* gr = new TGraphAsymmErrors(5, x, y, exl, exh, eyl, eyh);
  
  TLatex* cc_bs_label = new TLatex(gr->GetX()[0], gr->GetY()[0], " CC (B_{s})");
  TLatex* cc_bd_label = new TLatex(gr->GetX()[1], gr->GetY()[1], " CC (B_{d})");
  TLatex* cw_label = new TLatex(gr->GetX()[2], gr->GetY()[2], " CW");
  TLatex* ww_label = new TLatex(gr->GetX()[3], gr->GetY()[3], " WW");
  TLatex* cn_label = new TLatex(gr->GetX()[4], gr->GetY()[4], " CN");
  
  gr->GetListOfFunctions()->Add(cc_bs_label);
  gr->GetListOfFunctions()->Add(cc_bd_label);
  gr->GetListOfFunctions()->Add(cw_label);
  gr->GetListOfFunctions()->Add(ww_label);
  gr->GetListOfFunctions()->Add(cn_label);
  
  gr->SetTitle(title);
  gr->SetMarkerStyle(kOpenCircle);
  gr->SetMarkerColor(4);
  gr->Draw("AP");

  c1.Print(filename + ".eps");
  
  TFile accuracy_plot_file(filename + ".root", "RECREATE");
  gr->Write();
  accuracy_plot_file.Close();
  
  return;
}

void estimateSignalFitPerformance()
{
  //open the ROOT efficiency file
  TFile ROOTFile(efficiencyFile.c_str());
  if (!ROOTFile.IsOpen()) {
    cerr << "Error opening file " << efficiencyFile << ".\n";
    return;
  }

  //get numBackgroundFail, numBackgroundPass, numSignalAll, and efficiency
  RooFitResult* fitResult = NULL;
  ROOTFile.GetObject("PhotonToIDEB/unbinned/probe_eta_bin0__probe_nJets05_bin0__gaussPlusLinear/fitresults", fitResult);
  Double_t efficiency = 0.0;
  Double_t efficiencyError = 0.0;
  Double_t numBackgroundFail = 0.0;
  Double_t numBackgroundFailError = 0.0;
  Double_t numBackgroundPass = 0.0;
  Double_t numBackgroundPassError = 0.0;
  Double_t numSignalAll = 0.0;
  Double_t numSignalAllError = 0.0;
  if (fitResult != NULL) {
    RooRealVar* theEfficiency = (RooRealVar*)fitResult->floatParsFinal().find("efficiency");
    efficiency = theEfficiency->getVal();
    efficiencyError = theEfficiency->getError();
    RooRealVar* theNumBackgroundFail = (RooRealVar*)fitResult->floatParsFinal().find("numBackgroundFail");
    numBackgroundFail = theNumBackgroundFail->getVal();
    numBackgroundFailError = theNumBackgroundFail->getError();
    RooRealVar* theNumBackgroundPass = (RooRealVar*)fitResult->floatParsFinal().find("numBackgroundPass");
    numBackgroundPass = theNumBackgroundPass->getVal();
    numBackgroundPassError = theNumBackgroundPass->getError();
    RooRealVar* theNumSignalAll = (RooRealVar*)fitResult->floatParsFinal().find("numSignalAll");
    numSignalAll = theNumSignalAll->getVal();
    numSignalAllError = theNumSignalAll->getError();
  }
  else {
    cerr << "Error getting RooFitResult PhotonToIDEB/unbinned/probe_eta_bin0__probe_nJets05_bin0__gaussPlusLinear/fitresults from file ";
    cerr << efficiencyFile << ".\n";
  }

  //get integrals of tag-pass and tag-fail distributions
  TCanvas* fitCanvas = NULL;
  ROOTFile.GetObject("PhotonToIDEB/unbinned/probe_eta_bin0__probe_nJets05_bin0__gaussPlusLinear/fit_canvas", fitCanvas);
  Double_t tagPassIntegral = 0;
  Double_t tagFailIntegral = 0;
  if (fitCanvas != NULL) {
    fitCanvas->cd(1);
    RooHist* tagPassDistribution = NULL;
    tagPassDistribution = (RooHist*)((TCanvas*)fitCanvas->GetPrimitive("fit_canvas_1"))->GetPrimitive("h_data");
    fitCanvas->cd(2);
    RooHist* tagFailDistribution = NULL;
    tagFailDistribution = (RooHist*)((TCanvas*)fitCanvas->GetPrimitive("fit_canvas_2"))->GetPrimitive("h_data");
    RooHist* blah = NULL;
    blah = (RooHist*)((TCanvas*)fitCanvas->GetPrimitive("fit_canvas_3"))->GetPrimitive("h_data");
    cerr << blah->Integral() << endl;
    if ((tagPassDistribution != NULL) && (tagFailDistribution != NULL)) {
      tagPassIntegral = tagPassDistribution->Integral()*/*1.796*/1.844;
      tagFailIntegral = tagFailDistribution->Integral()*/*1.796*/1.844;
    }
    else cerr << "Error: could not get RooPlots.\n";
  }
  else {
    cerr << "Error getting TCanvas PhotonToIDEB/unbinned/probe_eta_bin0__probe_nJets05_bin0__gaussPlusLinear/fit_canvas from file ";
    cerr << efficiencyFile << ".\n";
  }

  //close file
  ROOTFile.Close();

  //subtract fitted background from integral
  Double_t tagPassNumBkgSubtractedEvts = tagPassIntegral - numBackgroundPass;
  Double_t tagPassNumBkgSubtractedEvtsError = numBackgroundPassError;
  Double_t tagFailNumBkgSubtractedEvts = tagFailIntegral - numBackgroundFail;
  Double_t tagFailNumBkgSubtractedEvtsError = numBackgroundFailError;

  //calculate fitted signal
  Double_t tagPassNumFittedSignal = efficiency*numSignalAll;
  Double_t tagPassNumFittedSignalError = tagPassNumFittedSignal*sqrt(((efficiencyError*efficiencyError)/(efficiency*efficiency)) + 
								     ((numSignalAllError*numSignalAllError)/(numSignalAll*numSignalAll)));
  Double_t tagFailNumFittedSignal = (1.0 - efficiency)*numSignalAll;
  Double_t tagFailNumFittedSignalError = tagFailNumFittedSignal*
    sqrt(((efficiencyError*efficiencyError)/((1.0 - efficiency)*(1.0 - efficiency))) + 
	 ((numSignalAllError*numSignalAllError)/(numSignalAll*numSignalAll)));

  //calculate difference between signal fit result and background subtracted integral
  Double_t tagPassDifference = tagPassNumBkgSubtractedEvts - tagPassNumFittedSignal;
  Double_t tagPassDifferenceError = sqrt(tagPassNumBkgSubtractedEvtsError*tagPassNumBkgSubtractedEvtsError + 
					 tagPassNumFittedSignalError*tagPassNumFittedSignalError);
  Double_t tagFailDifference = tagFailNumBkgSubtractedEvts - tagFailNumFittedSignal;
  Double_t tagFailDifferenceError = sqrt(tagFailNumBkgSubtractedEvtsError*tagFailNumBkgSubtractedEvtsError + 
					 tagFailNumFittedSignalError*tagFailNumFittedSignalError);

  //compare signal fit result to background subtracted integral
  cout << "Tag pass signal fit: " << tagPassNumFittedSignal << " +/- " << tagPassNumFittedSignalError << endl;
  cout << "Tag pass background subtracted integral: " <<  tagPassNumBkgSubtractedEvts << " +/- " << tagPassNumBkgSubtractedEvtsError;
  cout << endl;
  cout << "Difference: " << tagPassDifference << " +/- " << tagPassDifferenceError << endl;
  cout << "Tag fail signal fit: " << tagFailNumFittedSignal << " +/- " << tagFailNumFittedSignalError << endl;
  cout << "Tag fail background subtracted integral: " <<  tagFailNumBkgSubtractedEvts << " +/- " << tagFailNumBkgSubtractedEvtsError;
  cout << endl;
  cout << "Difference: " << tagFailDifference << " +/- " << tagFailDifferenceError << endl;
//.........这里部分代码省略.........

//.........这里部分代码省略.........

  cout << "Number of events in data sig: " << data->numEntries() << endl;
  cout << "Number of events in data bkg: " << data_bkg->numEntries() << endl;
  
  // Initialize tree to save toys to 
  TTree* results = new TTree("results","toy results");
  
  double fa3,fa3Error, fa3Pull;
  double sigFrac,sigFracError, sigFracPull;
  double significance;

  results->Branch("fa3",&fa3,"fa3/D");
  results->Branch("fa3Error",&fa3Error,"fa3Error/D");
  results->Branch("fa3Pull",&fa3Pull,"fa3Pull/D");
  results->Branch("sigFrac",&sigFrac,"sigFrac/D");
  results->Branch("sigFracError",&sigFracError,"sigFracError/D");
  results->Branch("sigFracPull",&sigFracPull,"sigFracPull/D");
  results->Branch("significance",&significance,"significance/D");

  //---------------------------------

  RooDataSet* toyData;
  RooDataSet* toyData_bkgOnly;
  int embedTracker=nEvts*counter;
  int embedTracker_bkg=TMath::Ceil(nEvts/3.75*counter);
  RooArgSet *tempEvent;

  RooFitResult *toyfitresults;
  RooFitResult *toyfitresults_sigBkg;
  RooFitResult *toyfitresults_bkgOnly;
  RooRealVar *r_fa3;
  RooRealVar *r_sigFrac;

  for(int i = 0 ; i<nToys ; i++){
    cout <<i<<"<-----------------------------"<<endl;
    //if(toyData) delete toyData;
    toyData = new RooDataSet("toyData","toyData",RooArgSet(*kd));
    toyData_bkgOnly = new RooDataSet("toyData_bkgOnly","toyData_bkgOnly",RooArgSet(*kd));

    if(nEvts+embedTracker > data->sumEntries()){
      cout << "Playground::generate() - ERROR!!! Playground::data does not have enough events to fill toy!!!!  bye :) " << endl;
      toyData = NULL;
      abort();
      return 0;
    }
    if(nEvts+embedTracker_bkg > data_bkg->sumEntries()){
      cout << "Playground::generate() - ERROR!!! Playground::data does not have enough events to fill toy!!!!  bye :) " << endl;
      toyData = NULL;
      abort();
      return 0;
    }

    for(int iEvent=0; iEvent<nEvts; iEvent++){
      if(iEvent==1)
	cout << "generating event: " << iEvent << " embedTracker: " << embedTracker << endl;
      tempEvent = (RooArgSet*) data->get(embedTracker);
      toyData->add(*tempEvent);
      embedTracker++;
    }
    if(bkg){
      for(int iEvent=0; iEvent<nEvts/3.75; iEvent++){
	if(iEvent==1)
	  cout << "generating bkg event: " << iEvent << " embedTracker bkg: " << embedTracker_bkg << endl;
	tempEvent = (RooArgSet*) data_bkg->get(embedTracker_bkg);
	toyData->add(*tempEvent);
	toyData_bkgOnly->add(*tempEvent);
	embedTracker_bkg++;
      }
    }

    if(bkg)
      toyfitresults =model.fitTo(*toyData,Save());
    else
      toyfitresults =modelSignal.fitTo(*toyData,Save());

    //cout<<toyfitresults<<endl;
    r_fa3 = (RooRealVar *) toyfitresults->floatParsFinal().find("fa3");

    fa3 = r_fa3->getVal();
    fa3Error = r_fa3->getError();
    fa3Pull = (r_fa3->getVal() - fa3Val) / r_fa3->getError();
    if(sigFloating){
      r_sigFrac = (RooRealVar *) toyfitresults->floatParsFinal().find("BoverTOT");
      sigFrac = 1-r_sigFrac->getVal();
      sigFracError = r_sigFrac->getError();
      sigFracPull = (1-r_sigFrac->getVal() - sigFracVal) / r_sigFrac->getError();
    }
    // fill TTree
    results->Fill();
  }

  char nEvtsString[100];
  sprintf(nEvtsString,"_%iEvts_%iiter",nEvts, counter);

  // write tree to output file (ouputFileName set at top)
  TFile *outputFile = new TFile("embeddedToys1DKD_fa3Corr_WithBackgDetEffects_"+sampleName[mySample]+nEvtsString+".root","RECREATE");
  results->Write();
  outputFile->Close();

}

prepDataFiles(){
//	TDirectory *theDr = (TDirectory*) myFile->Get("eleIDdir");///denom_pt/fit_eff_plots");
	//theDr->ls();
	int myIndex;	
	
	TSystemDirectory dir(thePath, thePath);
	TSystemFile *file;
	TString fname;
	TIter next(dir.GetListOfFiles());
	while ((file=(TSystemFile*)next())) {
		fname = file->GetName();
		if (fname.BeginsWith("TnP")&& fname.Contains("mc")) {
	
			ofstream myfile;

			TFile *myFile = new TFile(fname);
			TIter nextkey(myFile->GetListOfKeys());
			TKey *key;
			while (key = (TKey*)nextkey()) {
				TString theTypeClasse = key->GetClassName();
				TString theNomClasse = key->GetTitle();
				if ( theTypeClasse == "TDirectoryFile"){
					TDirectory *theDr = (TDirectory*) myFile->Get(theNomClasse);
					TIter nextkey2(theDr->GetListOfKeys());
					TKey *key2;
					while (key2 = (TKey*)nextkey2()) {
						TString theTypeClasse2 = key2->GetClassName();
						TString theNomClasse2 = key2->GetTitle();	
						myfile.open (theNomClasse2+".info");
						if ( theTypeClasse == "TDirectoryFile"){
							cout << "avant " << endl;
							TDirectory *theDr2 = (TDirectory*) myFile->Get(theNomClasse+"/"+theNomClasse2);
							cout << "apres " << endl;
							TIter nextkey3(theDr2->GetListOfKeys());
							TKey *key3;
							while (key3 = (TKey*)nextkey3()) {
								TString theTypeClasse3 = key3->GetClassName();
								TString theNomClasse3 = key3->GetTitle();	
								if ((theNomClasse3.Contains("FromMC"))) {

									TString localClasse3 = theNomClasse3;
									localClasse3.ReplaceAll("__","%");
									cout << "apres " << localClasse3 << endl;
									TObjArray* listBin = localClasse3.Tokenize('%');
									TString first = ((TObjString*)listBin->At(0))->GetString();
									TString second = ((TObjString*)listBin->At(2))->GetString();
									myfile << first;
									myfile << " " << second << " ";
									cout << "coucou la on va récupérer le rooFitResult " << endl;

									RooFitResult *theResults = (RooFitResult*) myFile->Get(theNomClasse+"/"+theNomClasse2+"/"+theNomClasse3+"/fitresults");
									theResults->Print();
									RooArgList theParam = theResults->floatParsFinal();
									int taille = theParam.getSize();
									for (int m = 0 ; m < taille ; m++){
										cout << "m=" << m << endl;
									RooAbsArg *theArg = (RooAbsArg*) theParam.at(m);
									RooAbsReal *theReal = (RooAbsReal*) theArg;
										myfile << theReal->getVal() << " " ;
									}		
															
									myfile << "\n";

								}
							}
						}
						myfile.close();

					}
			
				}
			}
			delete myFile;
		}
	
	}

}

//.........这里部分代码省略.........
  w->factory("SIMUL::McModel(McPassing,pass=McModelP,fail=McModelF)");
  // simultaneous fir for data
  w->factory("SIMUL::DataModel(DataPassing,pass=DataModelP,fail=DataModelF)");
  w->Print("V");
  w->saveSnapshot("clean", w->allVars());

  w->loadSnapshot("clean");

  /****************** sim fit to soup **************************/

  ///////////////////////////////////////////////////////////////
  TFile *f = new TFile("dummySoup.root","RECREATE");
  TTree* cutTreeSoupP = fullTreeSoup->CopyTree(Form("(%s>0.5 && %s && signalPFChargedHadrCands<1.5)",category_.c_str(),bin_.c_str()));
  TTree* cutTreeSoupF = fullTreeSoup->CopyTree(Form("(%s<0.5 && %s && signalPFChargedHadrCands<1.5)",category_.c_str(),bin_.c_str()));
 
  RooDataSet McDataP("McDataP","dataset pass for the soup", RooArgSet(*mass), Import( *cutTreeSoupP ) );
 
  RooDataSet McDataF("McDataF","dataset fail for the soup", RooArgSet(*mass), Import( *cutTreeSoupF ) );
 
  RooDataHist McCombData("McCombData","combined data for the soup", RooArgSet(*mass), Index(*(w->cat("McPassing"))), Import("pass", *(McDataP.createHistogram("histoP",*mass)) ), Import("fail",*(McDataF.createHistogram("histoF",*mass)) ) ) ;

  RooPlot* McFrameP    = 0;
  RooPlot* McFrameF    = 0;
  RooRealVar* McEffFit = 0;

  if(makeSoupFit_){

    cout << "**************** N bins in mass " << w->var("mass")->getBins() << endl;

    RooFitResult* ResMcCombinedFit = w->pdf("McModel")->fitTo(McCombData, Extended(1), Minos(1), Save(1),  SumW2Error( SumW2_ ), Range(xLow_,xHigh_), NumCPU(4) /*, ExternalConstraints( *(w->pdf("ConstrainMcNumBkgF")) )*/ );
    test->cd(Form("bin%f",binCenter_));
    ResMcCombinedFit->Write("McFitResults_Combined");

    RooArgSet McFitParam(ResMcCombinedFit->floatParsFinal());
    McEffFit     = (RooRealVar*)(&McFitParam["McEfficiency"]);
    RooRealVar* McNumSigFit  = (RooRealVar*)(&McFitParam["McNumSgn"]);
    RooRealVar* McNumBkgPFit = (RooRealVar*)(&McFitParam["McNumBkgP"]);
    RooRealVar* McNumBkgFFit = (RooRealVar*)(&McFitParam["McNumBkgF"]);

    McFrameP = mass->frame(Bins(24),Title("MC: passing sample"));
    McCombData.plotOn(McFrameP,Cut("McPassing==McPassing::pass"));
    w->pdf("McModel")->plotOn(McFrameP,Slice(*(w->cat("McPassing")),"pass"), ProjWData(*(w->cat("McPassing")),McCombData), LineColor(kBlue),Range(xLow_,xHigh_));
    w->pdf("McModel")->plotOn(McFrameP,Slice(*(w->cat("McPassing")),"pass"), ProjWData(*(w->cat("McPassing")),McCombData), Components("McSignalPdfP"), LineColor(kRed),Range(xLow_,xHigh_));
    w->pdf("McModel")->plotOn(McFrameP,Slice(*(w->cat("McPassing")),"pass"), ProjWData(*(w->cat("McPassing")),McCombData), Components("McBackgroundPdfP"), LineColor(kGreen),Range(xLow_,xHigh_));
    
    McFrameF = mass->frame(Bins(24),Title("MC: failing sample"));
    McCombData.plotOn(McFrameF,Cut("McPassing==McPassing::fail"));
    w->pdf("McModel")->plotOn(McFrameF,Slice(*(w->cat("McPassing")),"fail"), ProjWData(*(w->cat("McPassing")),McCombData), LineColor(kBlue),Range(xLow_,xHigh_));
    w->pdf("McModel")->plotOn(McFrameF,Slice(*(w->cat("McPassing")),"fail"), ProjWData(*(w->cat("McPassing")),McCombData), Components("McSignalPdfF"), LineColor(kRed),Range(xLow_,xHigh_)); 
    w->pdf("McModel")->plotOn(McFrameF,Slice(*(w->cat("McPassing")),"fail"), ProjWData(*(w->cat("McPassing")),McCombData), Components("McBackgroundPdfF"), LineColor(kGreen),Range(xLow_,xHigh_)); 
  }
  
  ///////////////////////////////////////////////////////////////

  /****************** sim fit to data **************************/

  ///////////////////////////////////////////////////////////////
  TFile *f2 = new TFile("dummyData.root","RECREATE");
  TTree* cutTreeDataP = fullTreeData->CopyTree(Form("(%s>0.5 && %s && signalPFChargedHadrCands<1.5)",category_.c_str(),bin_.c_str()));
  TTree* cutTreeDataF = fullTreeData->CopyTree(Form("(%s<0.5 && %s && signalPFChargedHadrCands<1.5)",category_.c_str(),bin_.c_str()));
 
  RooDataSet DataDataP("DataDataP","dataset pass for the soup", RooArgSet(*mass), Import( *cutTreeDataP ) );
  RooDataSet DataDataF("DataDataF","dataset fail for the soup", RooArgSet(*mass), Import( *cutTreeDataF ) );
  RooDataHist DataCombData("DataCombData","combined data for the soup", RooArgSet(*mass), Index(*(w->cat("DataPassing"))), Import("pass",*(DataDataP.createHistogram("histoDataP",*mass))),Import("fail",*(DataDataF.createHistogram("histoDataF",*mass)))) ;

  RooFitResult* ResDataCombinedFit = w->pdf("DataModel")->fitTo(DataCombData, Extended(1), Minos(1), Save(1),  SumW2Error( SumW2_ ), Range(xLow_,xHigh_), NumCPU(4));

//.........这里部分代码省略.........

        mZ = new RooFormulaVar("mZ", "TMath::Sqrt(2*@0+2*@1+2*@2+2*@3+2*@4+2*@[email protected]*@[email protected]*@7)", RooArgList(p1D2,p1Dph1,p2Dph1,p1Dph2,p2Dph2,ph1Dph2, m1, m2));
        RelBW = new RooGenericPdf("RelBW","1/( pow(mZ*mZ-bwMean*bwMean,2)+pow(mZ,4)*pow(bwGamma/bwMean,2) )", RooArgSet(*mZ,bwMean,bwGamma) );
//        PDFRelBW = new RooProdPdf("PDFRelBW", "PDFRelBW", RooArgList(gauss1_lep, gauss2_lep, gauss1_gamma, gauss2_gamma, *RelBW));

        }

     //true shape
     RooRealVar sg("sg", "sg", sgVal_);
     RooRealVar a("a", "a", aVal_);
     RooRealVar n("n", "n", nVal_);

     RooCBShape CB("CB","CB",*mZ,bwMean,sg,a,n);
     RooRealVar f("f","f", fVal_);

     RooRealVar mean("mean","mean",meanVal_);
     RooRealVar sigma("sigma","sigma",sigmaVal_);
     RooRealVar f1("f1","f1",f1Val_);

     RooAddPdf *RelBWxCB;
     RelBWxCB = new RooAddPdf("RelBWxCB","RelBWxCB", *RelBW, CB, f);
     RooGaussian *gauss;
     gauss = new RooGaussian("gauss","gauss",*mZ,mean,sigma);
     RooAddPdf *RelBWxCBxgauss;
     RelBWxCBxgauss = new RooAddPdf("RelBWxCBxgauss","RelBWxCBxgauss", *RelBWxCB, *gauss, f1);

     RooProdPdf *PDFRelBWxCBxgauss;
     PDFRelBWxCBxgauss = new RooProdPdf("PDFRelBWxCBxgauss","PDFRelBWxCBxgauss", 
                                     RooArgList(gauss1_lep, gauss2_lep, *RelBWxCBxgauss) );


    //make fit
    RooArgSet *rastmp;
    rastmp = new RooArgSet(pTRECO1_lep, pTRECO2_lep);
/*
    if(input.nFsr == 1) {
      rastmp = new RooArgSet(pTRECO1_lep, pTRECO2_lep, pTRECO1_gamma);
      }

    if(input.nFsr == 2) {
      rastmp = new RooArgSet(pTRECO1_lep, pTRECO2_lep, pTRECO1_gamma, pTRECO2_gamma);
      }
*/
    RooDataSet* pTs = new RooDataSet("pTs","pTs", *rastmp);
    pTs->add(*rastmp);

    RooFitResult* r;
    if (mass4lRECO_ > 140) {

       r = PDFRelBW->fitTo(*pTs,RooFit::Save(),RooFit::PrintLevel(-1));

       } else {

              r = PDFRelBWxCBxgauss->fitTo(*pTs,RooFit::Save(),RooFit::PrintLevel(-1));

              }
    //save fit result
    const TMatrixDSym& covMatrix = r->covarianceMatrix();
    const RooArgList& finalPars = r->floatParsFinal();

    for (int i=0 ; i<finalPars.getSize(); i++){
 
        TString name = TString(((RooRealVar*)finalPars.at(i))->GetName());
        if(debug_) cout<<"name list of RooRealVar for covariance matrix "<<name<<endl;

    }

    int size = covMatrix.GetNcols();
    output.covMatrixZ.ResizeTo(size,size);
    output.covMatrixZ = covMatrix;
    
    output.pT1_lep = pTMean1_lep.getVal();
    output.pT2_lep = pTMean2_lep.getVal();
    output.pTErr1_lep = pTMean1_lep.getError();
    output.pTErr2_lep = pTMean2_lep.getError();
/*
    if (input.nFsr >= 1) {

       output.pT1_gamma = pTMean1_gamma.getVal();
       output.pTErr1_gamma = pTMean1_gamma.getError();
    
       }

    if (input.nFsr == 2) {

       output.pT2_gamma = pTMean2_gamma.getVal();
       output.pTErr2_gamma = pTMean2_gamma.getError();

       }
*/
    delete rastmp;
    delete pTs;
    delete PDFRelBW;
    delete mZ;
    delete RelBW;
    delete RelBWxCB;
    delete gauss;
    delete RelBWxCBxgauss;
    delete PDFRelBWxCBxgauss;
}

//.........这里部分代码省略.........
//     MyLimit result(plInt->IsInInterval(exp_sig),
    MyLimit result(exp_sig_val>lowLim&&exp_sig_val<uppLim,lowLim,uppLim);
    // std::cout << "isIn " << result << std::endl;
    delete plInt;
//     delete modelConfig;
    return result;
  }

  // use FeldmaCousins (takes ~20 min)  
  if ( method == FeldmanCousinsMethod ) {
    FeldmanCousins fc(*data, modelConfig);
    fc.SetConfidenceLevel(0.95);
    //number counting: dataset always has 1 entry with N events observed
    fc.FluctuateNumDataEntries(false); 
    fc.UseAdaptiveSampling(true);
    fc.SetNBins(100);
    PointSetInterval* fcInt = NULL;
    fcInt = (PointSetInterval*) fc.GetInterval(); // fix cast
    double lowLim = fcInt->LowerLimit(*wspace->var("s"));
    double uppLim = fcInt->UpperLimit(*wspace->var("s"));
//     double exp_sig_val = wspace->var("s")->getVal();
    cout << "Feldman Cousins interval on s = [" << lowLim << " " << uppLim << endl;
    // std::cout << "isIn " << result << std::endl;
    MyLimit result(exp_sig_val>lowLim&&exp_sig_val<uppLim,
		   fcInt->LowerLimit(*wspace->var("s")),fcInt->UpperLimit(*wspace->var("s")));
    delete fcInt;
    return result;
  }


  // use BayesianCalculator (only 1-d parameter of interest, slow for this problem)  
  if ( method == BayesianMethod ) {
    BayesianCalculator bc(*data, modelConfig);
    bc.SetConfidenceLevel(0.95);
    bc.SetLeftSideTailFraction(0.5);
    SimpleInterval* bInt = NULL;
    if( wspace->set("poi")->getSize() == 1)   {
      bInt = bc.GetInterval();
      if ( draw ) {
	TCanvas* c = new TCanvas("Bayesian");
	// the plot takes a long time and print lots of error
	// using a scan it is better
	bc.SetScanOfPosterior(50);
	RooPlot* bplot = bc.GetPosteriorPlot();
	bplot->Draw();
      }
      cout << "Bayesian interval on s = [" << 
	bInt->LowerLimit( ) << ", " <<
	bInt->UpperLimit( ) << "]" << endl;
      // std::cout << "isIn " << result << std::endl;
      MyLimit result(bInt->IsInInterval(exp_sig),
		     bInt->LowerLimit(),bInt->UpperLimit());
      delete bInt;
      return result;
    } else {
    cout << "Bayesian Calc. only supports on parameter of interest" << endl;
    return MyLimit();
    }
  }


  // use MCMCCalculator  (takes about 1 min)
  // Want an efficient proposal function, so derive it from covariance
  // matrix of fit
  if ( method == MCMCMethod ) {
    RooFitResult* fit = wspace->pdf("model")->fitTo(*data,Save());
    ProposalHelper ph;
    ph.SetVariables((RooArgSet&)fit->floatParsFinal());
    ph.SetCovMatrix(fit->covarianceMatrix());
    ph.SetUpdateProposalParameters(kTRUE); // auto-create mean vars and add mappings
    ph.SetCacheSize(100);
    ProposalFunction* pf = ph.GetProposalFunction();
    
    MCMCCalculator mc(*data, modelConfig);
    mc.SetConfidenceLevel(0.95);
    mc.SetProposalFunction(*pf);
    mc.SetNumBurnInSteps(100); // first N steps to be ignored as burn-in
    mc.SetNumIters(100000);
    mc.SetLeftSideTailFraction(0.5); // make a central interval
    MCMCInterval* mcInt = NULL;
    mcInt = mc.GetInterval();
    MCMCIntervalPlot mcPlot(*mcInt); 
    mcPlot.Draw();
    cout << "MCMC interval on s = [" << 
      mcInt->LowerLimit(*wspace->var("s") ) << ", " <<
      mcInt->UpperLimit(*wspace->var("s") ) << "]" << endl;
    // std::cout << "isIn " << result << std::endl;
    MyLimit result(mcInt->IsInInterval(exp_sig),
		   mcInt->LowerLimit(*wspace->var("s")),mcInt->UpperLimit(*wspace->var("s")));
    delete mcInt;
    return result;
  }
  

  t.Print();

//   delete modelConfig;
  return MyLimit();

}

  void  accessParams(int n, double array[] ){
    /*
    RooArgList fl=fitresult->floatParsFinal();
    RooArgList cl=fitresult->constPars();
    fl.Print();
    cl.Print();
    int fn=fl.getSize();
    int cn=cl.getSize();
    */    //    mean,area
    RooRealVar* par;
    RooRealVar* pas;
    if (n==1){
                     par = (RooRealVar*) fitresult->floatParsFinal().find("mean1");
    if (par==NULL){  par = (RooRealVar*) fitresult->constPars().find("mean1"); }
                     pas = (RooRealVar*) fitresult->floatParsFinal().find("area1");
    if (pas==NULL){  pas = (RooRealVar*) fitresult->constPars().find("area1"); }
    //    printf( "   mean1 == %f +- %f\n", par->getVal() ,  par->getError()  );
    array[0]= par->getVal();     array[1]= par->getError();
    array[2]= pas->getVal();     array[3]= pas->getError();
    }
    if (npeaks<2) return;

    if (n==2){
                     par = (RooRealVar*) fitresult->floatParsFinal().find("mean2");
    if (par==NULL){  par = (RooRealVar*) fitresult->constPars().find("mean2"); }
                     pas = (RooRealVar*) fitresult->floatParsFinal().find("area2");
    if (pas==NULL){  pas = (RooRealVar*) fitresult->constPars().find("area2"); }
    //    printf( "   mean2 == %f +- %f\n", par->getVal() ,  par->getError()  );
    array[0]= par->getVal();     array[1]= par->getError();
     array[2]= pas->getVal();     array[3]= pas->getError();
   }
    if (npeaks<3) return;
    
    if (n==3){
                     par = (RooRealVar*) fitresult->floatParsFinal().find("mean3");
    if (par==NULL){  par = (RooRealVar*) fitresult->constPars().find("mean3"); }
                     pas = (RooRealVar*) fitresult->floatParsFinal().find("area3");
    if (pas==NULL){  pas = (RooRealVar*) fitresult->constPars().find("area3"); }
    //    printf( "   mean2 == %f +- %f\n", par->getVal() ,  par->getError()  );
    array[0]= par->getVal();     array[1]= par->getError();
    array[2]= pas->getVal();     array[3]= pas->getError();
    }
    if (npeaks<4) return;
    
    if (n==4){
                     par = (RooRealVar*) fitresult->floatParsFinal().find("mean4");
    if (par==NULL){  par = (RooRealVar*) fitresult->constPars().find("mean4"); }
                     pas = (RooRealVar*) fitresult->floatParsFinal().find("area4");
    if (pas==NULL){  pas = (RooRealVar*) fitresult->constPars().find("area4"); }
    //    printf( "   mean2 == %f +- %f\n", par->getVal() ,  par->getError()  );
    array[0]= par->getVal();     array[1]= par->getError();
     array[2]= pas->getVal();     array[3]= pas->getError();
   }
    if (npeaks<5) return;
    
    if (n==5){
                     par = (RooRealVar*) fitresult->floatParsFinal().find("mean5");
    if (par==NULL){  par = (RooRealVar*) fitresult->constPars().find("mean5"); }
                     pas = (RooRealVar*) fitresult->floatParsFinal().find("area5");
    if (pas==NULL){  pas = (RooRealVar*) fitresult->constPars().find("area5"); }
    //    printf( "   mean2 == %f +- %f\n", par->getVal() ,  par->getError()  );
    array[0]= par->getVal();     array[1]= par->getError();
     array[2]= pas->getVal();     array[3]= pas->getError();
   }
    if (npeaks<6) return;
    


    return;
  }//----------access