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

void voigtian(RooDataSet *dataset, RooDataSet *dataset2, RooRealVar &variable, RooPlot *fitFrame, RooBinning b, double rangeMin, double rangeMax, vector <double> &fitParameters)
{
        RooRealVar mean("mean","mean",0.0,-0.1,0.1);
        RooRealVar sigma("sigma","sigma",0.5,0.0,1.0);
        RooRealVar width("width","width",0.5,0.0,1.0);

        RooVoigtian * pdf = new RooVoigtian("pdf","Voigtian",variable,mean,sigma,width);

        dataset->plotOn(fitFrame,Name("myhist"),Binning(b),DataError(RooAbsData::SumW2));

        RooFitResult * res = pdf->fitTo(*dataset, Range(rangeMin, rangeMax),Save(),SumW2Error(kTRUE));
	res->Print();

	//minNll = res->minNll();
        
	pdf->plotOn(fitFrame,Name("mycurve"));


	fitParameters.push_back(3); // nb of fit parameters

	fitParameters.push_back(mean.getVal());
	fitParameters.push_back(mean.getError());

	fitParameters.push_back(sigma.getVal());
        fitParameters.push_back(sigma.getError());
	
	fitParameters.push_back(width.getVal());
        fitParameters.push_back(width.getError());		
	

}

void rf314_paramfitrange()
{

  // D e f i n e   o b s e r v a b l e s   a n d   d e c a y   p d f 
  // ---------------------------------------------------------------

  // Declare observables
  RooRealVar t("t","t",0,5) ;
  RooRealVar tmin("tmin","tmin",0,0,5) ;

  // Make parameterized range in t : [tmin,5]
  t.setRange(tmin,RooConst(t.getMax())) ;

  // Make pdf
  RooRealVar tau("tau","tau",-1.54,-10,-0.1) ;
  RooExponential model("model","model",t,tau) ;



  // C r e a t e   i n p u t   d a t a 
  // ------------------------------------

  // Generate complete dataset without acceptance cuts (for reference)
  RooDataSet* dall = model.generate(t,10000) ;

  // Generate a (fake) prototype dataset for acceptance limit values
  RooDataSet* tmp = RooGaussian("gmin","gmin",tmin,RooConst(0),RooConst(0.5)).generate(tmin,5000) ;

  // Generate dataset with t values that observe (t>tmin)
  RooDataSet* dacc = model.generate(t,ProtoData(*tmp)) ;



  // F i t   p d f   t o   d a t a   i n   a c c e p t a n c e   r e g i o n
  // -----------------------------------------------------------------------

  RooFitResult* r = model.fitTo(*dacc,Save()) ;


 
  // P l o t   f i t t e d   p d f   o n   f u l l   a n d   a c c e p t e d   d a t a 
  // ---------------------------------------------------------------------------------

  // Make plot frame, add datasets and overlay model
  RooPlot* frame = t.frame(Title("Fit to data with per-event acceptance")) ;
  dall->plotOn(frame,MarkerColor(kRed),LineColor(kRed)) ;
  model.plotOn(frame) ;
  dacc->plotOn(frame) ;

  // Print fit results to demonstrate absence of bias
  r->Print("v") ;


  new TCanvas("rf314_paramranges","rf314_paramranges",600,600) ;
  gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.6) ; frame->Draw() ;

  return ;
}

//_________________________________________________
void TestJeffreysGaussSigma(){
  // this one is VERY sensitive
  // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
  //   and you get really bizzare shapes
  // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
  //   and the PDF falls off too fast at high sigma
  RooWorkspace w("w");
  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])");
  w.factory("n[100,.1,2000]");
  w.factory("ExtendPdf::p(g,n)");
  //  w.var("sigma")->setConstant();
  w.var("mu")->setConstant();
  w.var("n")->setConstant();
  w.var("x")->setBins(301);

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;


  //  w.defineSet("poi","mu,sigma");
  //w.defineSet("poi","mu,sigma,n");
  w.defineSet("poi","sigma");
  w.defineSet("obs","x");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);

  const RooArgSet* temp = w.set("poi");
  pi.getParameters(*temp)->Print();
  //  return;

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","sqrt(2.)/sigma",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("sigma")->frame();
  pi.plotOn(plot);
  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
  plot->Draw();


}

  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

void JeffreysPriorDemo(){
  RooWorkspace w("w");
  w.factory("Uniform::u(x[0,1])");
  w.factory("mu[100,1,200]");
  w.factory("ExtendPdf::p(u,mu)");

  //  w.factory("Poisson::pois(n[0,inf],mu)");

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());
  //  RooDataHist* asimov2 = w.pdf("pois")->generateBinned(*w.var("n"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;

  w.defineSet("poi","mu");
  w.defineSet("obs","x");
  //  w.defineSet("obs2","n");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  //  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",10);

  //  JeffreysPrior pi2("jeffreys2","jeffreys",*w.pdf("pois"),*w.set("poi"),*w.set("obs2"));

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","1./sqrt(mu)",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("mu")->frame();
  //  pi.plotOn(plot, Normalization(1,RooAbsReal::Raw),Precision(.1));
  pi.plotOn(plot);
  //  pi2.plotOn(plot,LineColor(kGreen),LineStyle(kDotted));
  test->plotOn(plot,LineColor(kRed));
  plot->Draw();

}

//_________________________________________________
void TestJeffreysGaussMean(){
  RooWorkspace w("w");
  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,0,10])");
  w.factory("n[10,.1,200]");
  w.factory("ExtendPdf::p(g,n)");
  w.var("sigma")->setConstant();
  w.var("n")->setConstant();

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;

  //  w.defineSet("poi","mu,sigma");
  w.defineSet("poi","mu");
  w.defineSet("obs","x");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  //  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);

  const RooArgSet* temp = w.set("poi");
  pi.getParameters(*temp)->Print();

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","1",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("mu")->frame();
  pi.plotOn(plot);
  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
  plot->Draw();


}

RooFitResult* fitter_Zee(TH1D *hist){

  RooRealVar  Zmassvar("Zmassvar","Zmassvar", 82, 100);
  RooDataHist *datahist = new RooDataHist("data","Z Mass",Zmassvar,hist);
  RooPlot *Zmassvarframe = Zmassvar.frame(Name("Zmassvarframe"),Title(hist->GetTitle())) ;
  datahist->plotOn(Zmassvarframe);
  
  RooRealVar alpha  ("alpha"  ,        "alpha" , 0.005,0.001,0.1); 
  RooRealVar n      ("n"      ,            "n" , 1,0.001,10); 
  RooRealVar cbmean ("cbmean" ,       "cbmean" , 1, 0.8, 1.2);
  RooRealVar cbsigma("cbsigma",      "cbsigma" , 0.01, 0.001, 0.2);

  RooRealVar bwmean("bwmean","bwmean",91,85,95);
  RooRealVar bwsigma("bwsigma","bwsigma",3,2,4);
  RooRealVar expoconst("expoconst","expoconst",-0.1,-0.5,0);

  RooCBShape cball  ("cball"  , "crystal ball" , Zmassvar, cbmean, cbsigma, alpha, n);
  RooBreitWigner bw("bw","breit wigner",Zmassvar,bwmean,bwsigma);

  RooFFTConvPdf cballXbw("cballXbw","cball (X) bw",Zmassvar,bw,cball);
  RooExponential expo("expo", "exponential", Zmassvar, expoconst);

  RooRealVar frac("frac","frac",0.1,0.001,0.2);

  RooAddPdf Zshapemodel("Zshapemodel","expo + cball (X) bw",RooArgList(expo,cballXbw),frac);


  RooFitResult *fitres =Zshapemodel.fitTo(*datahist,Range(82,100),Save());   
  Zshapemodel.plotOn(Zmassvarframe,LineColor(kBlue));

  Zmassvarframe->Draw();

  fitres->Print();

  return fitres;

};

///
/// Test PDF implementation.
/// Performs a fit to the minimum.
///
bool PDF_Abs::test()
{
	bool quiet = false;
	if(quiet) RooMsgService::instance().setGlobalKillBelow(ERROR);
	fixParameters(observables);
	floatParameters(parameters);
	setLimit(parameters, "free");
	RooFormulaVar ll("ll", "ll", "-2*log(@0)", RooArgSet(*pdf));
	RooMinuit m(ll);
	if(quiet) m.setPrintLevel(-2);
	m.setNoWarn();
	m.setLogFile("/dev/zero");
	m.setErrorLevel(1.0);
	m.setStrategy(2);
	// m.setProfile(1);
	m.migrad();
	RooFitResult *f = m.save();
	bool status = !(f->edm()<1 && f->status()==0);
	if(!quiet) f->Print("v");
	delete f;
	if(quiet) RooMsgService::instance().setGlobalKillBelow(INFO);
	if(!quiet) cout << "pdf->getVal() = " << pdf->getVal() << endl;
	return status;
}

void SHyFT::fit(bool verbose)
{
  if ( verbose ) cout << "Fitting" << endl;
  RooArgSet nllset;
  for(unsigned int i=0;i<bins_.size();++i) {
    nllset.add(*bins_[i]->nll());
  }
  
  RooAddition nllsum("nllsum","nllsum",nllset);
  RooMinuit m(nllsum);

  if ( verbose ) m.setVerbose(kTRUE);
  else {
    m.setVerbose(kFALSE);
    m.setPrintLevel(-1);
  }
  m.migrad();
  m.hesse();
  //  m.minos();
  if ( verbose ) {
    RooFitResult * f = m.save();
    f->Print("v");
  }
}

//.........这里部分代码省略.........
  RooFormulaVar* nSigPass = new RooFormulaVar("nSigPass", "eff*numSignal", RooArgList(*eff,*numSignal));
  RooFormulaVar* nSigFail = new RooFormulaVar("nSigFail", "(1.0-eff)*numSignal", RooArgList(*eff,*numSignal));

  RooRealVar* nBkgPass = new RooRealVar("nBkgPass","nBkgPass", 100, 0.0, 10000000);
  if(SetPassBkgZero){
    nBkgPass = new RooRealVar("nBkgPass","nBkgPass", 0., 0., 0.);
  }


  RooRealVar* nBkgFail = new RooRealVar("nBkgFail","nBkgFail", 100, 0.0, 10000000);
  if(SetFailBkgZero){
    nBkgFail = new RooRealVar("nBkgFail","nBkgFail", 0.,0.,0.);
  }



  RooArgList componentsPass(*signalShapePdfPass,*bkgShapePdf);
  RooArgList yieldsPass(*nSigPass, *nBkgPass);

  RooArgList componentsFail(*signalShapePdfFail,*bkgShapePdff);
  // RooArgList componentsFail(BWxCB,*bkgShapePdff);


  RooArgList yieldsFail(*nSigFail, *nBkgFail);


   RooAddPdf pdfPass("pdfPass","extended sum pdf", componentsPass, yieldsPass);
   RooAddPdf pdfFail("pdfFail","extended sum pdf", componentsFail, yieldsFail);


   // The total simultaneous fit ...
   RooSimultaneous totalPdf("totalPdf","totalPdf", sample);
   totalPdf.addPdf(pdfPass,"Pass");
   totalPdf.Print();
   totalPdf.addPdf(pdfFail,"Fail");
   totalPdf.Print();


   ifstream readinfail(datafile_fail,ios::in);
   float mm;
   int ndataFailPeak = 0;
   while(readinfail.good()){
    if( readinfail.eof()) break;
    readinfail>>mm;

    if( mm> 80 && mm <100){
      ndataFailPeak ++;
    }

   }


   // ********* Do the Actual Fit ********** //

   RooFitResult *fitResult = totalPdf.fitTo(*data,RooFit::Save(true),
					    RooFit::Extended(true), RooFit::PrintLevel(-1));

   fitResult->Print("v");


  double numerator = nSigPass->getVal();
  double nfails    = nSigFail->getVal();
  double denominator = numerator + nfails;

  cout<<"num/den: "<< numerator <<" "<< denominator <<endl;

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

//.........这里部分代码省略.........
  RooHistPdf* signalShapeFailTemplatePdf = new RooHistPdf("signalShapeFailTemplatePdf", "signalShapeFailTemplatePdf", Mass, *dataHistFail, 1);

  RooFFTConvPdf* signalShapePassPdf = new RooFFTConvPdf("signalShapePassPdf","signalShapePassPdf"  , Mass, *signalShapePassTemplatePdf,*PassSignalResolutionFunction,2);
  RooFFTConvPdf* signalShapeFailPdf = new RooFFTConvPdf("signalShapeFailPdf","signalShapeFailPdf"  , Mass, *signalShapeFailTemplatePdf,*FailSignalResolutionFunction,2);


  // Now define some efficiency/yield variables  
  RooFormulaVar* NumSignalPass = new RooFormulaVar("NumSignalPass", "ParEfficiency*ParNumSignal", RooArgList(*ParEfficiency,*ParNumSignal));
  RooFormulaVar* NumSignalFail = new RooFormulaVar("NumSignalFail", "(1.0-ParEfficiency)*ParNumSignal", RooArgList(*ParEfficiency,*ParNumSignal));


  //*********************************************************************************************
  //
  // Create Background PDFs
  //
  //*********************************************************************************************
  RooExponential* bkgPassPdf = new RooExponential("bkgPassPdf","bkgPassPdf",Mass, *ParPassBackgroundExpCoefficient);
  RooExponential* bkgFailPdf = new RooExponential("bkgFailPdf","bkgFailPdf",Mass, *ParFailBackgroundExpCoefficient);


 //*********************************************************************************************
  //
  // Create Total PDFs
  //
  //*********************************************************************************************
  RooAddPdf pdfPass("pdfPass","pdfPass",RooArgList(*signalShapePassPdf,*bkgPassPdf), RooArgList(*NumSignalPass,*ParNumBkgPass));
  RooAddPdf pdfFail("pdfFail","pdfFail",RooArgList(*signalShapeFailPdf,*bkgFailPdf), RooArgList(*NumSignalFail,*ParNumBkgFail));

  // PDF for simultaneous fit
   RooSimultaneous totalPdf("totalPdf","totalPdf", sample);
   totalPdf.addPdf(pdfPass,"Pass");
//    totalPdf.Print();
   totalPdf.addPdf(pdfFail,"Fail");
   totalPdf.Print();


   //*********************************************************************************************
  //
  // Perform Fit
  //
  //*********************************************************************************************
   RooFitResult *fitResult = 0;

  // ********* Fix with Migrad first ********** //  
   fitResult = totalPdf.fitTo(*dataCombined, RooFit::Save(true), 
                              RooFit::Extended(true), RooFit::PrintLevel(-1));
   fitResult->Print("v");


//   // ********* Fit With Minos ********** //  
//    fitResult = totalPdf.fitTo(*dataCombined, RooFit::Save(true), 
//                               RooFit::Extended(true), RooFit::PrintLevel(-1), RooFit::Minos());
//    fitResult->Print("v");




//   // ********* Fix Mass Shift and Fit For Resolution ********** //  
//    ParPassSignalMassShift->setConstant(kTRUE); 
//    ParFailSignalMassShift->setConstant(kTRUE); 
//    ParPassSignalResolution->setConstant(kFALSE); 
//    ParFailSignalResolution->setConstant(kFALSE); 
//    fitResult = totalPdf.fitTo(*dataCombined, RooFit::Save(true), 
//    RooFit::Extended(true), RooFit::PrintLevel(-1));
//    fitResult->Print("v");

//.........这里部分代码省略.........
      histomaximum = histo[i][k]->GetBinContent(histo[i][k]->GetMaximumBin());
      
      // Declare observable x
      RooRealVar x("x","x",eneMIN,eneMAX) ;
      eneMIN = histo[i][k]->GetMaximumBin() - RMS; eneMAX = histo[i][k]->GetMaximumBin() + RMS;
      // std::cout << histo[i][k]->GetMaximumBin() << std::endl;
      // eneMIN = 50; eneMAX = 200.;
      x.setRange("fitRegion1",eneMIN  ,  eneMAX );

      RooDataHist dh("dh","dh",x,Import(*histo[i][k])) ;
      
      RooPlot* frame = x.frame(Title(filename[i]), Bins(1000)) ;
      dh.plotOn(frame,MarkerColor(2),MarkerSize(0.9),MarkerStyle(21));  //this will show histogram data points on canvas 
      dh.statOn(frame);  //this will display hist stat on canvas
      
      RooRealVar mean("mean","mean",histo_mean, eneMIN, eneMAX);
      RooRealVar width("width","width",RMS, eneMIN, eneMAX);
      RooRealVar sigma("sigma","sigma",RMS, eneMIN, eneMAX);
      
       RooGaussian gauss("gauss","gauss",x,mean,sigma);
      // RooExponential expo("expo", "exponential PDF", x, lambda);
      
      // Construct landau(t,ml,sl) ;
      //RooLandau landau("lx","lx",x,mean,sigma) ;
  
      //RooBreitWigner BW1("BW1","Breit Wigner theory",x,mean,sigma);
      RooVoigtian bwandgaus("bwandgaus", "Breit Wigner convoluted with gauss",x,mean,width,sigma);
      //RooCBShape cryBall1("cryBall1","Crystal Ball resolution model", x, CBmean, CBsigma, CBalpha, CBn) ;
      //RooFFTConvPdf bwxCryBall1("bwxCryBall1", "FFT Conv CryBall and BW", x, BW1, cryBall1); 
      // RooBreitWigner gauss("gauss","gauss",x,mean,sigma);
      // RooVoigtian gauss("gauss","gauss",x,mean,width,sigma);

      // RooFitResult* filters = gauss.fitTo(dh, NumCPU(4), Range(("fitRegion"+IntToString(i)).c_str()), Save(true), SumW2Error(kTRUE), PrintLevel(-1));
       RooFitResult* filters = gauss.fitTo(dh, NumCPU(4), Range("fitRegion1"), Save(true), SumW2Error(kTRUE), PrintLevel(-1));
      //RooFitResult* filters = bwxCryBall1.fitTo(dh, NumCPU(4), Range("fitRegion1"), Save(true), SumW2Error(kTRUE), PrintLevel(-1));
      // RooFitResult* filters = BW1.fitTo(dh, NumCPU(4), Range("fitRegion1"), Save(true), SumW2Error(kTRUE), PrintLevel(-1));
      //RooFitResult* filters = bwandgaus.fitTo(dh, NumCPU(4), Range("fitRegion1"), Save(true), SumW2Error(kTRUE), PrintLevel(-1));
      // RooFitResult* filters = sum.fitTo(dh, NumCPU(4), Range("fitRegion1"), Save(true), SumW2Error(kTRUE), PrintLevel(-1));
      // RooFitResult* filters = gauss.fitTo(dh, NumCPU(4), Save(true), SumW2Error(kTRUE), PrintLevel(-1));

      //RooFitResult* filters = landau.fitTo(dh, NumCPU(4), Range("fitRegion1"), Save(true), SumW2Error(kTRUE), PrintLevel(-1));
      gauss.plotOn(frame,LineColor(4));//this will show fit overlay on canvas 
      gauss.paramOn(frame); //this will display the fit parameters on canvas
      
      filters->Print();
      
      // // Access list of final fit parameter values
      // std::cout << "final value of floating parameters" << std::endl ;
      // filters->floatParsFinal().Print("s") ;
      
      //mpc = fit->GetParameter(1) - mpshift * fit->GetParameter(2); 
      mpc = mean.getVal();

      maxv.push_back(mpc);
      maxvr.push_back( (mean.getError()) );
      eneerr.push_back(0.);

      thesigma.push_back( (sigma.getVal()) );
      thesigmaerr.push_back( (sigma.getError()) );
   
      // std::cout << "File " << filename[i] << " Eneer " << k << " MAXV "<< mpc << std::endl;

      // can[i][k]->Print(canname + ".root",".root");
      // can[i][k]->Close();
  
  	// // Draw all frames on a canvas

void CreateTemplatesForW(TString CAT, TString CUT)
{
  gROOT->ForceStyle();
  
  RooMsgService::instance().setSilentMode(kTRUE);
  RooMsgService::instance().setStreamStatus(0,kFALSE);
  RooMsgService::instance().setStreamStatus(1,kFALSE);

  TFile *infMC   = TFile::Open("Histo_TT_TuneCUETP8M1_13TeV-powheg-pythia8.root");
  TFile *infData = TFile::Open("Histo_JetHT.root");
    
  TH1F *hMC,*hData;
  RooDataHist *roohMC,*roohData;
  RooAddPdf *signal,*qcd;
  
  TString VAR,TAG;
  float XMIN,XMAX;

  VAR = "mW";
  TAG = CUT+"_"+CAT;
  XMIN = 20.;
  XMAX = 160.;

  RooWorkspace *w = new RooWorkspace("w","workspace");

  //---- define observable ------------------------
  RooRealVar *x = new RooRealVar("mW","mW",XMIN,XMAX);
  w->import(*x);
  //---- first do the data template ---------------
  
  hData = (TH1F*)infData->Get("boosted/h_mW_"+CUT+"_0btag");
  hData->Rebin(5);
  roohData = new RooDataHist("roohistData","roohistData",RooArgList(*x),hData);  

  //---- QCD -----------------------------------
  RooRealVar bBkg0("qcd_b0","qcd_b0",0.5,0,1);
  RooRealVar bBkg1("qcd_b1","qcd_b1",0.5,0,1);
  RooRealVar bBkg2("qcd_b2","qcd_b2",0.5,0,1);
  RooRealVar bBkg3("qcd_b3","qcd_b3",0.5,0,1);

  RooBernstein qcd1("qcd_brn","qcd_brn",*x,RooArgList(bBkg0,bBkg1,bBkg2,bBkg3));

  RooRealVar mQCD("qcd_mean2" ,"qcd_mean2",40,0,100);
  RooRealVar sQCD("qcd_sigma2","qcd_sigma2",20,0,50);
  RooGaussian qcd2("qcd_gaus" ,"qcd_gaus",*x,mQCD,sQCD);

  RooRealVar fqcd("qcd_f1","qcd_f1",0.5,0,1);

  qcd = new RooAddPdf("qcd_pdf","qcd_pdf",qcd1,qcd2,fqcd);

  //---- plots ---------------------------------------------------
  TCanvas *canB = new TCanvas("Template_W_QCD_"+CUT+"_"+CAT,"Template_W_QCD_"+CUT+"_"+CAT,900,600);

  RooFitResult *res = qcd->fitTo(*roohData,RooFit::Save());
  res->Print();
  RooPlot *frameB = x->frame();
  roohData->plotOn(frameB);
  qcd->plotOn(frameB);
  qcd->plotOn(frameB,RooFit::Components("qcd_brn"),RooFit::LineColor(kRed),RooFit::LineWidth(2),RooFit::LineStyle(2));
  qcd->plotOn(frameB,RooFit::Components("qcd_gaus"),RooFit::LineColor(kGreen+1),RooFit::LineWidth(2),RooFit::LineStyle(2));
  frameB->GetXaxis()->SetTitle("m_{W} (GeV)");
  frameB->Draw();
  gPad->Update();
  canB->Print("plots/"+TString(canB->GetName())+".pdf");

  RooArgSet *parsQCD = (RooArgSet*)qcd->getParameters(roohData);
  parsQCD->setAttribAll("Constant",true);

  w->import(*qcd);

  //---- then do the signal templates -------------
  hMC = (TH1F*)infMC->Get("boosted/hWt_"+VAR+"_"+TAG);
  roohMC = new RooDataHist("roohistTT","roohistTT",RooArgList(*x),hMC);
  normMC = ((TH1F*)infMC->Get("eventCounter/GenEventWeight"))->GetSumOfWeights();
  hMC->Rebin(2);
      
  RooRealVar m("ttbar_mean","ttbar_mean",90,70,100);
  RooRealVar s("ttbar_sigma","ttbar_sigma",20,0,50);
    
  RooGaussian sig_core("ttbar_core","ttbar_core",*x,m,s);
  
  RooRealVar bSig0("ttbar_b0","ttbar_b0",0.5,0,1);
  RooRealVar bSig1("ttbar_b1","ttbar_b1",0.5,0,1);
  RooRealVar bSig2("ttbar_b2","ttbar_b2",0.5,0,1); 
  RooRealVar bSig3("ttbar_b3","ttbar_b3",0.5,0,1);
  RooRealVar bSig4("ttbar_b4","ttbar_b4",0.5,0,1);
  RooRealVar bSig5("ttbar_b5","ttbar_b5",0.5,0,1); 
  RooRealVar bSig6("ttbar_b6","ttbar_b6",0.5,0,1);
  RooRealVar bSig7("ttbar_b7","ttbar_b7",0.5,0,1);
  RooRealVar bSig8("ttbar_b8","ttbar_b8",0.5,0,1);

  RooBernstein sig_tail("ttbar_tail","ttbar_tail",*x,RooArgList(bSig0,bSig1,bSig2,bSig3,bSig4,bSig5,bSig6,bSig7,bSig8)); 

  RooRealVar fcore("ttbar_fcore","ttbar_fcore",0.5,0,1);

  signal = new RooAddPdf("ttbar_pdf","ttbar_pdf",RooArgList(sig_core,sig_tail),RooArgList(fcore));

  TCanvas *canS = new TCanvas("Template_W_TT_"+CAT+"_"+CUT,"Template_W_TT_"+CAT+"_"+CUT,900,600);

  res = signal->fitTo(*roohMC,RooFit::Save());
//.........这里部分代码省略.........

   void ws_constrained_profile3D( const char* wsfile = "rootfiles/ws-data-unblind.root",
                                   const char* new_poi_name = "n_M234_H4_3b",
                                   int npoiPoints = 20,
                                   double poiMinVal = 0.,
                                   double poiMaxVal = 20.,
                                   double constraintWidth = 1.5,
                                   double ymax = 10.,
                                   int verbLevel=0 ) {


     gStyle->SetOptStat(0) ;

     //--- make output directory.

     char command[10000] ;
     sprintf( command, "basename %s", wsfile ) ;
     TString wsfilenopath = gSystem->GetFromPipe( command ) ;
     wsfilenopath.ReplaceAll(".root","") ;
     char outputdirstr[1000] ;
     sprintf( outputdirstr, "outputfiles/scans-%s", wsfilenopath.Data() ) ;
     TString outputdir( outputdirstr ) ;


     printf("\n\n Creating output directory: %s\n\n", outputdir.Data() ) ;
     sprintf(command, "mkdir -p %s", outputdir.Data() ) ;
     gSystem->Exec( command ) ;


     //--- Tell RooFit to shut up about anything less important than an ERROR.
      RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR) ;



       if ( verbLevel > 0 ) { printf("\n\n Verbose level : %d\n\n", verbLevel) ; }


       TFile* wstf = new TFile( wsfile ) ;

       RooWorkspace* ws = dynamic_cast<RooWorkspace*>( wstf->Get("ws") );

       if ( verbLevel > 0 ) { ws->Print() ; }






       RooDataSet* rds = (RooDataSet*) ws->obj( "ra2b_observed_rds" ) ;

       if ( verbLevel > 0 ) {
          printf("\n\n\n  ===== RooDataSet ====================\n\n") ;
          rds->Print() ;
          rds->printMultiline(cout, 1, kTRUE, "") ;
       }





       ModelConfig* modelConfig = (ModelConfig*) ws->obj( "SbModel" ) ;
       RooAbsPdf* likelihood = modelConfig->GetPdf() ;

       RooRealVar* rrv_mu_susy_all0lep = ws->var("mu_susy_all0lep") ;
       if ( rrv_mu_susy_all0lep == 0x0 ) {
          printf("\n\n\n *** can't find mu_susy_all0lep in workspace.  Quitting.\n\n\n") ;
          return ;
       }





       //-- do BG only.
       rrv_mu_susy_all0lep->setVal(0.) ;
       rrv_mu_susy_all0lep->setConstant( kTRUE ) ;










       //-- do a prefit.

       printf("\n\n\n ====== Pre fit with unmodified nll var.\n\n") ;

       RooFitResult* dataFitResultSusyFixed = likelihood->fitTo(*rds, Save(true),Hesse(false),Minos(false),Strategy(1),PrintLevel(verbLevel));
       int dataSusyFixedFitCovQual = dataFitResultSusyFixed->covQual() ;
       if ( dataSusyFixedFitCovQual < 2 ) { printf("\n\n\n *** Failed fit!  Cov qual %d.  Quitting.\n\n", dataSusyFixedFitCovQual ) ; return ; }
       double dataFitSusyFixedNll = dataFitResultSusyFixed->minNll() ;

       if ( verbLevel > 0 ) {
          dataFitResultSusyFixed->Print("v") ;
       }

       printf("\n\n Nll value, from fit result : %.3f\n\n", dataFitSusyFixedNll ) ;

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